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Abstract 


Space— Time  Fluid— Structure  Interaction  Computation  of 
Flapping- Wing  Aerodynamics 

by 

Nikolay  M.  Rostov 


The  focus  of  this  thesis  is  computational  fluid-structure  interaction  (FSI)  analysis 
of  flapping-wing  aerodynamics  of  a  micro  aerial  vehicle  (MAV).  The  wing  motion 
and  deformation  data,  whether  prescribed  fully  or  partially,  is  from  an  actual  locust, 
extracted  from  high-speed,  multi-camera  video  recordings  of  the  locust  in  a  wind 
tunnel.  The  core  computational  FSI  technology  is  based  on  the  Deforming-Spatial- 
Domain/Stabilized  Space-Time  (DSD/SST)  formulation.  This  is  supplemented  with 
using  NURBS  basis  functions  in  temporal  representation  of  the  wing  and  mesh  mo¬ 
tion,  and  in  remeshing.  Here  we  use  the  version  of  the  DSD/SST  formulation  derived 
in  conjunction  with  the  variational  multiscale  (VMS)  method,  and  this  version  is 
called  ’’DSD/SST-VMST.”  The  structural  mechanics  computations  are  based  on  the 
Kir chhoff- Love  shell  model.  We  use  a  sequential  coupling  technique,  which  is  ap¬ 
plicable  to  some  classes  of  FSI  problem,  especially  those  with  temporally-periodic 
behavior.  We  show  that  sequential  coupling  performs  well  in  FSI  computation  of 
the  flapping-wing  aerodynamics  we  consider  here.  We  analyze  cases  where  the  MAV 
body  has  rolling,  pitching,  or  rolling  and  pitching  motion.  We  study  how  all  these 
influence  the  lift  and  thrust  generated  by  the  MAV. 
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Chapter  1 


Introduction 


A  micro  aerial  vehicle  (MAV)  is  an  aerial  vehicle  with  a  scale  generally  smaller  than 
60  cm.  There  are  many  current  and  potential  applications  for  MAVs:  exploration  of 
large  areas  with  minimal  human  monitoring,  mapping  of  geographic  areas,  creation  of 
3D  surface  models  based  on  image  capturing,  surveillance,  monitoring  of  hazardous 
situations,  reduction  of  human  exposure  to  toxic  materials,  monitoring  warehouse  in¬ 
ventory,  and  many  others.  All  of  these  rely  on  self-sufficient  devices  that  can  sustain 
long  periods  of  flight  without  a  need  to  recharge  or  refuel.  This  requires  an  energy- 
efficient  flight  system  and  it  is  only  natural  to  look  at  biological  examples  of  high 
efficiency  flight  at  small  scale.  Fixed-wing  aerodynamics  is  not  very  efficient  at  a 
small  scale  as  the  forward  velocity  is  not  enough  to  generate  the  aerodynamic  forces 
necessary  without  increasing  wing  area  too  much.  This  requires  having  wings  that 
move  faster  than  the  forward  velocity.  This  can  be  achieved  by  either  rotating  blades 
or  flapping  wings.  Flapping-wing  insects  are  known  to  be  very  energy  efficient.  Lo¬ 
custs,  in  particular,  can  fly  very  long  distances  over  several  days  on  a  small  amount  of 
calories.  Bio-inspired  flapping  wings  allow  a  high  energy  efficiency  and  maneuverabil¬ 
ity.  However,  to  be  able  to  understand  and  utilize  these  biological  flapping  patterns 
we  need  accurate  simulation  and  modeling  techniques.  Computational  modeling  of 
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flapping-wing  aerodynamics  can  allow  a  cost-effective  way  to  develop  future  MAVs. 
In  this  thesis  we  will  present  the  numerical  methods  used  for  MAV  design  including 
the  fluid-structure  interactions  (FSI)  involved  in  the  motion  and  deformation  of  the 
wings  in  flight. 

1.1  Project  description 

The  ability  to  accurately  predict  the  aerodynamic  forces  during  flapping  flight  allows 
flexibility  in  the  design  of  MAVs.  With  little  additional  cost  many  different  wing 
conhgurations  can  be  evaluated  with  no  physical  models  being  built.  In  addition, 
wing  interactions  can  be  studied  to  maximize  benehcial  aerodynamic  effects.  A  variety 
of  flight  maneuvers  can  be  studied,  including  different  angles  of  attack,  rolling,  and 
pitching.  The  rest  of  the  material  in  this  section  is  from  [2]. 

Space-time  (ST)  computational  analysis  of  flapping-wing  aerodynamics  of  an  ac¬ 
tual  locust  was  hrst  presented  in  [3] .  The  wing  motion  and  deformation  data  was  ex¬ 
tracted  from  high-speed,  multi-camera  video  recordings  of  the  locust  in  a  wind  tunnel, 
which  was  part  of  a  parallel,  experimental  research  carried  out  by  Professor  Fabrizio 
Gabbiani  and  Dr.  Raymond  Chan  at  Baylor  College  of  Medicine  in  Houston.  The 
ST  computational  analysis  was  based  on  the  Deforming-Spatial-Domain/Stabilized 
ST  (DSD/SST)  formulation  [4,  5,  6,  7,  8,  9,  10,  11],  which  is  the  core  computational 
method,  and  special  ST  techniques  targeting  flapping-wing  aerodynamics. 

The  DSD/SST  formulation  is  an  interface-tracking  (moving-mesh)  method  for 
computation  of  flows  with  moving  interfaces,  including  FSI.  The  stabilization  com¬ 
ponents  of  the  formulation  are  the  Streamline-Upwind/Petrov-Galerkin  (SUPG)  [12] 
and  Pressure-Stabilizing/Petrov-Galerkin  (PSPG)  [4,  13]  methods.  The  ST-VMS 
method  [10]  is  the  the  variational  multiscale  version  of  the  DSD/SST  method,  which 
was  originally  called  “DSD/SST-VMST”  (i.e.  the  version  with  the  VMS  turbulence 
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model)  in  [9] .  The  VMS  components  are  from  the  residual-based  VMS  method  given 
in  [14,  15,  16,  17].  The  original  DSD/SST  formulation  was  named  “DSD/SST-SUPS” 
in  [9]  (i.e.  the  version  with  the  SUPG/PSPG  stabilization),  which  was  also  called 
“ST-SUPS”  in  [11], 

In  interface-tracking  methods,  as  the  interface  moves  and  the  spatial  domain  oc¬ 
cupied  by  the  fluid  changes  its  shape,  the  mesh  moves  to  accommodate  this  shape 
change  and  to  follow  (i.e.  “track”)  the  interface.  The  Arbitrary  Lagrangian-Eulerian 
(ALE)  hnite  element  formulation  [18]  is  the  most  widely  used  moving-mesh  technique, 
with  increased  emphasis  on  ESI  in  recent  years  (see,  for  example,  [19,  20,  21,  22,  23, 
24,  25,  26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,  42,  43,  44,  45,  46, 
11,  47,  48,  49,  50,  51]). 

Moving-mesh  methods  require  mesh  update  methods.  Mesh  update  typically  con¬ 
sists  of  moving  the  mesh  for  as  long  as  possible  and  remeshing  as  needed.  With  the  key 
objectives  being  to  maintain  the  element  quality  near  solid  surfaces  and  to  minimize 
frequency  of  remeshing,  a  number  of  advanced  mesh  update  methods  [52,  53,  54,  8] 
were  developed  to  be  used  with  the  DSD/SST  method,  including  those  that  minimize 
the  deformation  of  the  layers  of  small  elements  placed  near  solid  surfaces. 

Though  not  as  popular  as  the  ALE  formulation,  the  DSD/SST  method  has  been 
applied  to  some  of  the  most  challenging  flow  problems  with  moving  interfaces.  The 
classes  of  problems  solved  include  the  free-surface  and  multi-fluid  flows  [4,  6,  52,  55, 
54],  fluid-object  interactions  [4,  5,  6,  55],  flows  with  solid  surfaces  in  fast,  linear  or 
rotational  relative  motion  [55,  54,  56,  57,  37,  58],  compressible  flows  [55],  fluid-particle 
interactions  [55,  54],  and  ESI  [59,  60,  61,  62,  8,  63,  64,  65,  66,  67,  68,  69,  70,  71,  72, 
9,  73,  74,  75,  76,  10,  77,  78,  79]. 

One  of  the  desirable  features  of  the  DSD/SST  formulation  is  that  for  the  temporal 
representation  in  ST  computations  we  can  use  higher-order  basis  functions  such  as 
NURBS,  which  have  been  used  very  effectively  as  spatial  basis  functions  in  fluid  and 


22 


4 


structural  mechanics  and  FSI  computations  [80,  21,  25,  81].  Using  NURBS  basis 
functions  for  the  temporal  representation  provides  a  more  accurate  representation  of 
the  surface  motion  and  deformation  and  a  more  robust  and  effective  way  of  moving 
the  mesh  and  remeshing.  This  is  at  the  core  of  the  special  ST  techniques  targeting 
flapping-wing  aerodynamics.  The  special  ST  techniques  were  demonstrated  in  [3]  with 
the  preliminary  computation  of  the  flapping-wing  aerodynamics  of  an  actual  locust.  A 
number  of  advances  in  the  special  ST  techniques  were  presented  in  [1].  The  advances 
include  a  more  accurate  and  realistic  representation  of  the  wing  geometry,  a  smoother 
temporal  representation  of  the  wing  motion  and  deformation,  a  temporally-periodic 
representation  of  the  wing  motion  and  deformation,  a  more  effective  way  of  generating 
and  moving  an  inner  mesh  around  the  locust,  and  a  better  starting  condition  for 
the  mesh  motion.  A  detailed  computational  analysis  of  bio-inspired  flapping-wing 
aerodynamics  of  an  MAV  was  presented  in  [82],  where  the  wing  shapes,  motions  and 
deformations  were  essentially  the  same  as  those  of  the  actual  locust  in  [1].  The  MAV 
body  was  a  simplified  and  streamlined  version  of  the  locust  body  in  [1] .  The  detailed 
computational  analysis  included  temporal  and  spatial  refinement  studies,  and  also 
using  different  combinations  of  wing  conhgurations  and  investigating  the  beneficial 
and  disruptive  interactions  between  the  wings  and  the  role  of  wing  camber  and  twist. 
A  condensed  version  of  the  material  from  [1],  concentrating  on  the  flapping-motion 
modeling  and  computations,  and  a  temporal-order  study  from  [82]  were  presented 
in  [83]. 

In  an  FSI  computation  with  a  moving-mesh  method,  the  FSI  coupling  technique 
determines  how  the  coupling  between  the  equation  blocks  representing  the  fluid  me¬ 
chanics,  structural  mechanics,  and  mesh  moving  equations  is  handled.  The  coupling 
techniques  used  in  ST  FSI  computations  (i.e.  FSI  computations  with  the  DST/SST 
method)  evolved  over  the  years  from  block-iterative  FSI  coupling  [84]  (see  [62,  8]  for 
the  terminology)  to  a  more  robust  version  of  block-iterative  coupling  [84,  62,  85]  and 
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to  quasi-direct  coupling  [62]  and  direct  coupling  [62]  techniques.  The  quasi-direct 
and  direct  coupling  techniques,  which  are  applicable  to  cases  with  nonmatching  fluid 
and  structure  meshes  at  the  interface,  yield  more  robust  algorithms  for  FSI  com¬ 
putations  where  the  structure  is  light,  such  as  parachute  FSI  computations.  New 
versions  of  the  quasi-direct  and  direct  coupling  techniques  with  upgraded  and  addi¬ 
tional  interface  projection  methods  [8,  64,  69,  72,  86]  have  a  substantially  increased 
robustness  in  ST  FSI  computations,  and  rendered  the  earlier  ST  FSI  solvers  obso¬ 
lete.  These  new  quasi-direct  and  direct  coupling  techniques  automatically  reduce  to 
“monolithic”  coupling  when  the  interface  has  matching  fluid  and  structure  meshes. 
Allowing  nonmatching  meshes  at  the  interface  substantially  increases  the  scope  of 
the  FSI  solver,  leading  to  success  in  FSI  modeling  of  challenging  problems,  such  as 
ringsail  spacecraft  parachutes  (see  [64,  65,  72,  86,  87,  73,  75,  77,  78]). 

The  sequentially-coupled  FSI  (SCFSI)  technique  was  introduced  in  [88,  89]  in 
the  context  of  arterial  FSI  computations.  It  is  applicable  to  some  classes  of  FSI 
problems,  especially  those  with  temporally-periodic  FSI  dynamics,  where  the  FSI 
coupling  can  be  achieved  by  a  sequence  of  alternating  fluid  and  structural  mechanics 
computations,  with  each  computation  taking  place  over  a  cycle  or  more  of  the  FSI 
dynamics.  Multiscale  versions  of  the  SCFSI  technique  increased  the  effectiveness  of 
ST  FSI  computations  in  arterial  fluid  mechanics  [68]  and  parachute  modeling  [72,  87]. 

In  this  paper  we  present  a  SCFSI  analysis  of  flapping-wing  aerodynamics  of  an 
MAV.  The  wing  motion  and  deformation  data,  whether  prescribed  fully  or  partially,  is 
the  same  as  those  of  the  actual  locust  in  [1].  The  body  is  the  same  as  the  one  in  [82]. 
In  the  fluid  mechanics  computations,  we  use  the  ST-VMS  method  in  combination 
with  the  ST-SUPS  method.  The  structural  mechanics  computations  are  mostly  based 
on  the  Kirchhoff-Love  shell  model  [90,  91,  33].  We  also  carry  out,  for  comparison 
purpose,  some  test  computations  with  the  membrane  model.  In  addition  to  the 
straight-flight  case,  we  analyze  cases  where  the  MAV  body  has  rolling,  pitching,  or 
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rolling  and  pitching  motion,  where  we  assume  that  the  MAV  somehow  performs  such 
maneuvers.  We  study  how  these  maneuvers  influence  the  lift  and  thrust. 

1.2  Overview 

In  Chapter  2  we  review  the  governing  equations  for  the  fluid  and  structural  mechanics 
models.  MAV  flight  is  at  small  Mach  numbers,  so  the  Navier-Stokes  equations  of 
incompressible  flows  are  used.  The  DSD/SST  and  DSD/SST-VMST  formulations 
are  reviewed. 

In  Chapter  3  we  review  the  special  techniques  used  for  the  representation  of  the 
wing  motion,  mesh  generation,  and  mesh  motion.  We  use  NURBS  basis  functions  in 
time,  which  allows  more  accurate  wing  motion  representation.  In  addition,  we  review 
the  techniques  for  enforcing  a  no-slip  boundary  condition  and  for  starting  the  fluid 
mechanics  computation  when  the  motion  is  described  using  NURBS  in  time. 

In  Chapter  4  we  describe  the  steps  involved  in  fluid  mechanics  computation  of  a 
flapping-wing  MAV  while  using  actual  locust  data. 

In  Chapter  5  we  describe  the  mesh  generation  for  the  fluid  mechanics  computation 
of  the  MAV  with  locust  wings.  The  aerodynamic  forces  generated  are  calculated  and 
reported. 

In  Chapter  6  we  describe  the  SCFSI  technique  used  for  MAV  computations.  For 
the  computations  in  this  chapter  we  use  the  Kirchhoff-Love  shell  model  for  the  struc¬ 
ture.  The  wing  deformation  pattern  is  based  on  prescribing  the  leading  edge  of  the 
wing  with  the  motion  from  the  locust  data  and  solving  for  the  rest  of  the  wing  defor¬ 
mation.  Lift  and  thrust  results  are  presented  and  compared  to  the  MAV  with  wing 
motion  coming  entirely  from  the  locust  data. 

In  Chapter  7  we  report  the  results  from  the  SCFSI  computation  of  the  MAV 
with  the  membrane  model  for  the  wings.  The  results  for  the  lift  and  thrust  for  the 
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membrane  case  are  compared  to  the  shell  case. 

In  Chapter  8  we  discuss  a  number  of  cases  where  we  prescribe  the  motion  of 
the  MAV.  First,  we  discuss  a  case  where  we  use  the  actual  locust  body  motion  and 
evaluate  the  effect  of  that  motion  on  the  aerodynamic  forces  generated.  Next,  we 
prescribe  large  rolling,  pitching,  and  combined  rolling  and  pitching  motions.  The  lift 
and  thrust  generated  in  these  cases  are  discussed.  The  SCFSI  technique  with  the 
shell  model  is  used  for  these  large  prescribed- mot  ion  cases. 

In  Chapter  9  we  evaluate  variations  in  the  flapping  motion.  We  alter  the  synchro¬ 
nization  between  the  forewings  and  hindwings  and  study  the  effect  on  the  aerody¬ 
namic  forces  generated.  We  also  look  at  the  effect  of  asymmetric  flapping,  where  the 
motions  of  the  left  and  right  wings  are  not  the  same.  The  SCFSI  technique  with  the 
shell  model  is  used  for  these  cases. 

In  Chapter  10  we  provide  our  concluding  remarks. 
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Chapter  2 

Governing  equations  and  finite 
element  formulations 

Flapping  flight  is  inherently  an  FSI  problem.  The  wing  motion  generates  aerody¬ 
namic  forces,  which  in  turn  deform  the  wings.  This  interaction  yields  the  complex 
deformations  and  fluid  flow  that  determine  the  dynamics  of  flapping  flight.  The  fluid 
mechanics  is  governed  by  the  Navier-Stokes  equations  of  incompressible  flows.  These 
equations  are  reviewed  in  Section  2.1.  The  governing  equations  of  motion  for  the 
structural  mechanics  are  presented  in  Section  2.2.  The  DSD/SST  formulation  is  re¬ 
viewed  in  2.3.  The  DSD/SST-VMST  formulation  is  described  in  Section  2.4.  The 
new-generation  ST  formulations,  emphasizing  temporal  interpolation  with  NURBS 
basis  functions,  are  described  in  Section  2.5.  The  semi-discrete  formulation  of  struc¬ 
tural  mechanics  is  described  in  Section  2.6.  The  material  in  this  chapter  is  from  [10]. 

2.1  Fluid  mechanics  governing  equations 

Let  Qt  C  be  the  spatial  domain  with  boundary  F*  at  time  t  G  (0,T).  The  sub¬ 
script  t  indicates  the  time-dependence  of  the  domain.  The  Navier-Stokes  equations 
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of  incompressible  flows  are  written  on  and  Vf  G  (0,  T)  as 

pl^^  +  u-Vu-f^  -  V  -a  =  0,  (2.1) 

V-u  =  0,  (2.2) 

where  p,  u  and  f  are  the  density,  velocity  and  the  external  force,  respectively.  The 
stress  tensor  cr  is  dehned  as  cr{p,  u)  =  —pl  +  2p,e{u),  with  e(u)  =  ((Vu)  +  (Vu)^)  /2. 
Here  p  is  the  pressure,  I  is  the  identity  tensor,  p,  =  pu  is  the  viscosity,  z/  is  the 
kinematic  viscosity,  and  e(u)  is  the  strain-rate  tensor.  The  essential  and  natural 
boundary  conditions  for  Eq.  (2.1)  are  represented  as  u  =  g  on  (rt)g  and  n  ■  <t  =  h  on 
(rt)h,  where  (Tjg  and  (rt)h  are  complementary  subsets  of  the  boundary  T^,  n  is  the 
unit  normal  vector,  and  g  and  h  are  given  functions.  A  divergence-free  velocity  held 
uo(x)  is  specihed  as  the  initial  condition. 


2.2  Structural  mechanics  governing  equations 


Let  hlf  C  be  the  spatial  domain  with  boundary  T^,  where  Uxd  =  2  for  shells 

and  membranes.  The  superscript  “s”  indicates  the  structure.  The  parts  of  T^  corre¬ 
sponding  to  the  essential  and  natural  boundary  conditions  are  represented  by  (rf)g 
and  (rf)j^.  The  equations  of  motion  are  written  as 


P 


<Py  ^  dy  ,, 


V  .  ff*  =  0, 


(2.3) 


where  p®,  y,  p,  and  cr®  are  the  material  density,  structural  displacement,  damping 
coefficient,  external  force  and  the  Cauchy  stress  tensor,  respectively.  The  stresses  are 
expressed  in  terms  of  the  second  Piola-Kirchoff  stress  tensor  S,  which  is  related  to 
the  Cauchy  stress  tensor  through  a  kinematic  transformation. 

For  the  happing-wing  MAV,  what  makes  one  structural  element  model  different 
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from  the  other  is  the  manner  in  which  S  is  defined.  These  definitions  can  be  fonnd 
in  [11], 

2.3  DSD/SST  formulation  of  fluid  mechanics 

A  ST  variational  formnlation  of  incompressible  flows  (see  for  example  [4,  5,  6,  7, 
8,  9])  is  written  over  a  seqnence  of  N  ST  slabs  where  Qn  is  the  slice  of  the 
ST  domain  between  the  time  levels  and  tn+i,  and  is  the  lateral  bonndary  of 
Qn  (see  Figure  2.1).  We  denote  the  trial  and  test  functions  spaces  for  the  velocity 


Figure  2.1:  ST  slab  in  an  abstract  representation. 

and  pressure  as  u  G  5u,  p  G  Sp,  w  G  Vu  and  q  G  Vp.  The  notation  (■)“  and 
(■)+  denotes  the  function  values  at  as  approached  from  below  and  above.  At 

each  time  step,  the  integrations  are  performed  over  Qn-  The  essential  and  natural 
boundary  conditions  are  enforced  over  {PQg  and  {PQh,  the  complementary  subsets  of 
the  lateral  boundary  of  the  ST  slab.  In  the  DSD/SST  method  [4,  5,  6,  7,  8,  9,  10,  11], 
the  ST  finite  element  interpolation  functions  are  continuous  within  a  ST  slab,  but 
discontinuous  from  one  ST  slab  to  another.  Each  Qn  is  decomposed  into  elements 
Q® ,  where  e  =  1,  2, . . . ,  {nei)n-  The  subscript  n  used  with  Uei  is  for  the  general  case 
where  the  number  of  ST  elements  may  change  from  one  ST  slab  to  another.  The 
hnite-dimensional  trial  and  test  functions  spaces  are  denoted  as 
and 

The  DSD/SST  formulation  (from  [7])  is  written  as  follows:  given  (u^)7,  find 
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e  (S^)n  and  e  (iS^  )„,  such  that  Vw^  G  (Vu)n  and  Vg^  G  (V^  )„: 


■  p 


'Qr, 


(9u^ 


+  •  Vu^  -  fM  dQ+  f  e(w^)  :  o-(/,  u^)  dQ 

/  ^  Qn 


/  w*  .  h*  dP  +  /  /V  ■u'‘dQ+  (w"):  .  p  ((u'-):  -  (u'-).:)  dSl 

^  (^n)h  ^ Qn  ^ ^n 

(«'ez)Ti 


^  (9w^ 


yx  /  “['^supgp  ( 

("^eOn  « 

yX  /  ^LsicV  ■  w'^pV  ■  u'"  dQ  =  0, 

e=l  •'Qn 


+  u"  ■  Vw"  +  TpsPoVj'"! .  lL(p",  u")  -  pCIdQ 


(2.4) 


where 


L(g\  w^)  =  p  -  V  ■  a(g^ 


w’^). 


(2.5) 


This  formulation  is  applied  to  all  ST  slabs  Qo,  Qi,  Q2,  •••,  Qn-i,  starting  with 
(u^)q  =  uq.  Here  Tsupg,  t'pspg  and  z/psic  are  the  SUPG,  PSPG  and  LSIG  stabilization 
parameters.  There  are  various  ways  of  dehning  these  parameters.  Here  we  provide 
the  dehnitions  given  in  [7]: 


'  SUPG 


'^SUGN12  — 


'Gugns  — 


+ 


Tc 


SUGN12 

ON, 


T 

I  o 


'^en 

E 

^a=l 
^RGN 

4z/  ’ 


“  +  ■  VNa 


dt 


(2.6) 

(2.7) 

(2.8) 


-1 


^RGN  —  2  (  ^  ^  |r  •  ^Na 


.  a=l 


Vllu'’ 


r  = 


Vllu^ 


'TspG  —  'GupG) 


(2.9) 

(2.10) 


and  in  [8]: 


_  II  ^  ^11^ 

^LSIC  =  'GuPG  II  U  V  II  , 


(2.11) 
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where  rien  is  the  number  of  (ST)  element  nodes  and  Na  is  the  ST  shape  function 
associated  with  the  ST  node  a.  As  an  alternative  to  the  construction  of  Tsupg  as 
given  by  Eqs.  (2.6)-(2.7),  another  option  was  introduced  in  [8].  In  that  option, 
t'supg  is  constructed  based  on  separate  dehnitions  for  the  advection-dominated  and 
transient-dominated  limits: 


T'supg 


'O'SUGNI 


'^SUGN2 


At 

T’ 


(2.12) 

(2.13) 

(2.14) 


where  is  the  mesh  velocity  and  At  is  the  time-step  size.  It  was  noted  in  [8]  that 
separating  rsuGNi2  into  its  advection-  and  transient-dominated  components  as  given 
by  Eqs.  (2.13)-(2.14)  is  equivalent  to  excluding  the  part  of  in  Eq.  (2.7), 

making  that  the  dehnition  for  Tsugni,  and  accounting  for  in  the  dehnition  for 

Augn2  given  by  Eq.  (2.14).  Here  ^  is  the  vector  of  element  coordinates.  For  more 
ways  of  calculating  Tsupg,  t'pspg  and  r'Lsic,  see  [92,  7,  93,  94,  95,  96,  97,  98,  99,  100, 
101,  102,  103,  104,  105,  106]. 


2.4  DSD/SST-VMST  (ST- VMS)  formulation  of  fluid 
mechanics 

The  ST  version  of  the  residual-based  VMS  (RBVMS)  method  [14,  15,  16,  17]  was 
introduced  in  [9]  and  compared  to  the  original  DSD/SST  formulation  [4,  5,  6,  7,  8]. 

We  repeat  that  derivation  and  comparison  in  this  section. 
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2.4.1  ST  variational  formulation 

In  deriving  the  variational  formulation,  we  use  the  version  of  Eq.  (2.1)  where  the 
advection  term  is  in  the  conservation- law  form.  We  multiply  that  version  of  Eq.  (2.1) 
and  Eq.  (2.2)  with  the  corresponding  test  functions  and  integrate  them  over  Qn'- 


(  w  ■  p  +'V  ■  (uu)  —  fj  dQ  —  (  w  ■  V  ■  trdQ  +  [  qV  ■  udQ  =  0.  (2.15) 

jQn  J  JQu  JQn 

We  integrate  by  parts  all  the  terms  except  for  the  external  force  and  enforce  the 
essential  and  natural  boundary  conditions.  That  gives  us  the  following  variational 
formulation;  End  u  G  iSu  and  p  G  Sp,  such  that  V  w  G  Vu  and  V  g  G  Vp 

[  ■  pu-^^dn  -  [  w+  ■  ^  ■  pudQ 

J  f!„+l  J  J  Qn 

—  (w  ■  pu)  (n  ■  v)  dP  +  (w  ■  pu)  (n  ■  u)  dP  —  /  Vw  ;  puudQ 

(-Pn)h  {^n)h  ^  Qn 

—  w  ■  pidQ  —  w  •  hdP  +  /  ^(w)  :  adQ  +  qn  -  udP 

^  Qn  ^  (Pn)h  ^  Qn  ^  Pn 


—  /  Vq  •  udQ  =  0. 


(2.16) 


'Qn 


2.4.2  Scale  separation 

In  the  variational  multiscale  techniques  [14,  15,  16,  17]  the  “coarse-scale”  and  “fine- 
scale”  are  separated  as  follows: 


‘5u  —  Sp  —  Sp  (B  Sp,  Vu  —  Vu  ©  Vu,  Vp  —  Vp  ©  Vp. 


(2.17) 
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The  coarse-scale  part  of  Eq.  (2.16)  is  written  as  follows: 


[  /  w+-pu^dn-[  ^-pudQ 

'0„+l  JUn  d  Qn 

-  (w  ■  pu)  {n-'v)dP+  (w  ■  pu)  {n  ■  u)  dP  —  /  Vw  ; 

^  (-Pn)h  ^  (-Pri)h  ^  Qn 


puudQ 


—  w  ■  pidQ  —  w  ■  hdP  +  /  ^(w)  ;  adQ  +  qn  -  udP 

'd  Qn  'd  (Pri)h  Qn  'd 

—  /  Vg  ■  ndQ  =  0. 

''  Qn 


(2.18) 


From  [14,  15,  16,  17],  the  hne-scale  solutions  are  represented  by  the  strong-form 
residuals  of  the  coarse-scale: 


u  = -^fm  (u,p) ,  p  = -pz/crc  (u) , 


(2.19) 


where 


fm  (u,p) 
rc  (u) 


p  +  u  ■  Vu  -  f  )  +Vp  -  2V  ■  pe(u). 


=  Vu, 


(2.20) 

(2.21) 


and  tm  and  uc  are  stabilization  parameters  closely  related  to  Tsupg/t'pspg  and  z/lsig- 

Remark  1  More  on  the  fine-seale  approximation  in  eonjunetion  with  the  Green’s 
operator  ean  he  found  in  [If,  15,  16,  17]. 


2.4.3  DSD/SST-VMST  (ST- VMS)  formulation 

In  this  subsection  we  write  from  [9]  the  derivation  of  the  DSD/SST-VMST  formulation 
(i.e.  the  version  with  the  variational  multiscale  turbulence  model). 
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Fine-scale  representation 

The  fine-scale  solutions  are  represented  over  each  element  from  Eq.  (2.19)  with  G 
(‘5u)„  and  e 

u'  =  -— rM^’/),  p'  =-pz/crcK).  (2.22) 

p 

Remark  2  When  the  polynomial  order  of  the  shape  funetions  is  less  than  two,  the 
last  term  in  Eq.  (2.20)  vanishes. 


There  are  various  ways  of  defining  tm  and  z/q-  For  tm  we  use  the  Tsupg  definition 
given  by  Eq.  (2.6).  For  vq,  we  consider  the  z^lsic  definition  given  by  Eq.  (2.11)  and 
the  definition  from  [17]: 


z/c  = 


^sd  o/-  ‘~\/- 

^  d^k  d^k 

2-^  ■ 

fc=l 


(2.23) 


We  evaluate  the  stabilization  parameters  at  ^  =  0. 


Coarse-scale  discretization 

Using  scale  separation  in  Eq.  (2.18)  for  velocity  and  pressure,  we  obtain: 


'Qn 


(w)„+i  ■  p  ((u)„+i  +  (u\+i)  dVt-  (w)+  ■  p  ((u)„  +  (u\)  dVt 


Ow  f 

— —  ■  p(u  -|-  u')dQ  +  (w  ■  p  (u  -|-  u'))  (n  ■  (u  -|-  u'  —  v))  dP 


dt 


/  Vw  :  p(u -|- u')(u -|- u')(i(5  —  /  w  ■  pfdQ  —  /  w -hdP 

^  Qn  ^  Qn  ^ 

/  £(w)  :  (a(p,  u) -|- ff')  dQ -l-  /  gn  ■  (u -|- u')dP  —  /  Vg  ■  (u -|- u')d(5  =  0. 

V  Qn  V  Pn  Qn 


(2.24) 


34 


16 


Here  a'  =  a  —  a  is  introduced  temporarily.  We  set  the  fine-scale  solution  to  zero 
at  the  spatial  and  temporal  boundaries,  use  the  assumption  e(w)  ;  2/iVu'  =  0  (see 
[107,  108]),  and  obtain  the  following  form: 


(w)_  1  ■  p(u)_  -  /  (w)+  ■  p(u)„  dn 


^n  +  1 


f  (9w  f 

/  —■p{u  +  u)dQ+  /  (w  ■  pu)  (n  ■  (u  -  v))  dP 

iQn  J(Pr.\ 


—  /  Vw  ;  p(u -|- u')(u -|- u')d(5  —  /  w  ■  pfdQ  —  /  w  ■  hdP 

^  Qn  ^  Qn  ^  (^n)h 

-|-  /  e(w)  :  ((r(p,  u)  —  p'l)  dQ -l-  /  gn  ■  udP  —  /  Vg  ■  (u -|- u')d(5  =  0.  (2.25) 

^  Qn  ^  Pn  ^  Qn 

We  collect  the  fine-scale  terms  in  one  place  and  write  the  global  integrations  as  sum 
of  element-level  integrals: 


f  Mn+l-  pHn+ld^-  [  ^  '  pudQ 

'0„  +  l  d  j  Qn 


/  (w  ■  pu)  (n  ■  (u  —  v))  dP  —  /  Vw  :  pu  udQ  —  /  w  ■  pfdQ 

(^n)h  ^  Qn  ^  Qn 


—  w  ■  hdP  +  /  £(w)  :  a{p,  u)dQ  +  qn  -  udP  —  /  Vg  ■  udQ 

^  {Pn)ii  ^  Qn  ^  Pn  ^  Qn 

inel)7 

-E 


e=l  ■'Qn 


(9w 


p-^  +  Vg  )  ■  u'  +  Vw  :  (p(u'u  +  uu'  +  u'u')  +  p'l) 


dQ  =  0.  (2.26) 
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Spatially  discretized  version  of  Eq.  (2.26)  is  written  as  follows:  find  G  and 

G  such  that  V  G  (V^)^  and  V  G  (Vp)^: 


■  p{u'^)^dn 


^n  +  l 


'Qn 


(9W^ 

~w 


pu^dQ 


+ 


f  (w^  ■  pu’^)  (n'^  ■  (u'^  -  v'^))  dP-  [  Vw^  ;  ■  pi^dQ 

^  {Pn)i^  ^  Qn  ^  Qn 

M  ihiTj  I  f  I  f  „hh  ^^hiry  f  T-7  h 


/  ■  h'^dP  +  /  £(w'^)  ;  o■(p^  u'^)dg  +  /  q^n'^  ■  u'^dP  -  Vq^  ■  u'^dQ 

^  {Pn)l^  Qn  ^  Pn  ^  Qn 

Ki)n  „  r.  ^  h  \ 

J  (p-^  +  VgM  ■  U  +  Vw'^  ;  (p(u'u'^  +  u'^u'  +  u'u')  +  p'l) 


e=l  “'Qn  L 


dQ  =  0. 


(2.27) 


Comparison  with  the  original  DSD/SST  formulation 

To  compare  DSD/SST- VMST  (ST-VMS)  to  a  slightly  modihed  DSD/SST-SUPS 
(where  the  advection  term  is  in  the  conservation-law  form),  we  further  rearrange 
the  terms  in  Eq.  (2.27): 


■  p 


'Qn 


du^- 


+  V  ■  (u'^u'^)  -  dQ+  [ 

/  ^  Qn 


£(w'')  :  u^)dQ 


w’^-h^dP+  g"V-u'^dg+/  (w'^)^p((u")+-(u");)dD 


>^Pn)h 


'Qn 


S  -^Qn 


P 


/ 


+  ■  Vw'*  +  Vq 


iP-el^n  p 

■  \idQ  —  /  V  ■  w^pdQ 

e=l  “'Qn 


pu'  ■  (Vw^)  ■  u'dg  =  0. 


(2.28) 


Remark  3  The  last  two  terms  correspond  to  the  cross  stress  and  Reynolds  stress. 


Remark  4  If  we  exclude  the  last  two  terms,  the  formulation  becomes  the  same  as  the 
modified  DSD/SST-SUPS  formulation  (where  the  advection  term  is  in  the  conservation- 
law  form)  under  the  conditions  Tpspg  =  Tsupg  and  vq  =  Vlsic-  The  6th  and  7th  terms 


are  the  SUPG/PSPG  and  LSIG  terms. 
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Equation  (2.28)  is  conservative  form  of  DSD/SST-VMST  method.  The  convective 
form  of  DSD/SST-VMST  is  directly  comparable  to  the  unmodihed  DSD/SST-SUPS. 
For  a  shortcut  derivation  of  the  convective  form,  we  integrate  by  parts  two  of  the 
pieces  in  the  5th  term  of  Eq.  (2.25): 


+ 


+ 


'Q 


'Q 


Vw  ;  puudQ  =  —  (w  ■  pu)  (n  ■  u)  dP  +  w  ■  p(u  ■  Vu)dQ 

^  (^ri)h  Qn 

(w  ■  pu)  V  ■  udQ,  (2.29) 

/  Vw  ;  pvLudQ  =  —  (w  ■  pu)  (n  ■  u')  dP  +  w  ■  p(u'  ■  Vu)dQ 

jQr,. 


/  (w  ■  pu)  V  ■  vidQ. 

^  Qn 


(2.30) 


The  first  term  on  the  right- hand- side  of  Eq.  (2.29)  cancels  the  same  term  with  opposite 
sign  in  Eq.  (2.25).  The  first  term  on  the  right-hand-side  of  Eq.  (2.30)  vanishes  because 
the  hne-scale  solution  is  zero  at  the  spatial  boundaries.  The  sum  of  the  last  terms 
of  Eqs.  (2.29)  and  Eq.  (2.30)  is  zero.  The  discrete  version  of  the  second  term  on  the 
right-hand-side  of  Eq.  (2.29)  replaces  in  the  first  term  of  Eq.  (2.28)  the  advection 
piece  in  the  conservation-law  form.  The  discrete  version  of  the  second  term  on  the 
right-hand-side  of  Eq.  (2.30)  replaces  the  8th  term  of  Eq.  (2.28).  With  all  that,  we 
obtain  the  convective  form  of  DSD/SST-VMST: 


■  p 


'Qn 


w 


du^' 


+  dQ+  [  e(w'^)  :  ff(p^  u’^)dQ 

/  Qn 


(nel)n 

-v 


'‘■h‘‘dP+  /  q'‘V-n>‘dQ+  /  (w*)*  .  p  ((u^y  -  (u*);)  rff! 


'Qn 


e=l  JQ^n  L 


\  dt 


+  u'"  ■  Vw'"  -h  Vg 


■  u'dQ  —  /  V  ■  w^p'dQ 

e=l  JQ'k 


{P^edn  p  (^ez)n  n 

+  J2  p^'' 

e=l  e=l 


(2.31) 
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2.5  New-generation  ST  formulations 

A  ST  basis  function  can  be  written  as  a  product  of  its  spatial  and  temporal  parts: 

N:  =  T“  (0)  iV,  (0  ,  a  =  1,  2, . . . ,  n,ns,  a  =  1,2,...,  n,nu  (2-32) 

where  6  G  [—1,  1]  is  the  temporal  element  coordinate,  and  rign  and  rient  are  the  number 
of  spatial  and  temporal  element  nodes  (see  Figure  2.2  for  examples  of  temporal  basis 
functions).  In  general,  the  values 

(p-  =  lim.  (j)’^  (t) ,  0+  =  lini  cp^  (t) ,  (2.33) 

do  not  need  to  be  equal  to  0”Ti  and  (pl^  (coefficients  of  the  basis  functions  and 


Figure  2.2:  Temporal  basis  functions. 

T^).  However,  for  the  cases  we  consider  here  the  basis  functions  are  interpolatory  at 
6  =  —1  and  0  =  1,  and  therefore  (p~  =  (p^Tl  and  <pp  =  <p\. 

As  pointed  out  in  [9],  different  components  (i.e.  unknowns),  and  the  corresponding 
test  functions,  can  be  discretized  with  different  sets  of  temporal  basis  functions.  This 
was  shown  in  [9]  by  introducing  a  secondary  mapping  0(^(6')  G  [—1,  1],  where  0^(6') 
is  a  strictly-increasing  function,  and  rewriting  the  generalized  ST  basis  function  for 
the  element  indices  (a,  a)  as 


(7V“)^  =  T“(0c(0))iV,(O. 


(2.34) 
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Here  C  indicates  the  component,  which  can  also  be  “t”.  We  can  also  use  different 
functions  T“  for  different  components,  which  is  in  fact  something  we  do  in  another 
context,  and  this  could  be  in  combination  with  the  secondary  mapping.  Prescribed 
and  unknown  variables  can  be  represented  over  different  ST  slabs,  because  we  only 
need  to  supply  the  prescribed  values  at  the  integration  points  of  the  ST  slab  used  in 
representing  the  unknown  variables.  Most  of  the  time  higher-order  basis  functions 
can  represent  complex  functions  with  fewer  number  of  control  points.  This  is  very 
helpful  in  decreasing  the  I/O  intensity,  such  as  in  computations  with  the  multiscale 
SCFSI  techniques  (see  Section  7  of  [10]). 


2.5.1  Mesh  representation 

In  general,  for  the  moving-mesh  and  FSI  cases,  the  temporal  basis  functions  for  the 
mesh  represent  the  mesh  path.  As  pointed  out  in  [9],  we  need  the  flexibility  to  specify 
the  speed  along  the  path.  Consider  the  mesh  position  vector  given  by 

x(0)  =  T“(0.(0))iVX,  (2.35) 


where  the  repeated  indices  imply  summation  over  the  applicable  range.  It  was  pro¬ 
posed  in  [9]  to  use  0x(6')  to  specify  the  velocity.  We  seek  a  function  0x  such  that 


- /V  X 

dt  de  d0x  “  “ 


(2.36) 


represents  the  desired  mesh  velocity  along  the  path.  This  approach  gives  us  the  option 
to  form  a  ST  parametric  space  beyond  the  NURBS  capability. 


39 


21 


2.6  Semi-discrete  formulation  of  structural  mechan¬ 
ics 


With  and  coming  from  appropriately  defined  trial  and  test  function  spaces, 
respectively,  the  semi-discrete  hnite  element  formulation  of  the  structural  mechanics 
equations  (see  [109,  110,  111])  is  written  as 


■P° 


dt"^ 


dVt  + 


■VP 


dt 


+  p^r)  dn. 


dVL  + 


6E'^  dQ  = 


(2.37) 


The  fluid  mechanics  forces  acting  on  the  structure  are  represented  by  vector  The 
above  formulation  is  for  structures  represented  by  a  shell  or  membrane  model.  The 
left-hand-side  terms  of  Eq.  (2.37)  are  referred  to  in  the  original  conhguration  and  the 
right-hand-side  terms  in  the  deformed  configuration  at  time  t.  Time  discretization  of 
Eq.  (2.37)  is  based  on  the  Generalized-a  method  [112]. 
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Chapter  3 


Special  techniques 


The  representation  of  the  motion  of  the  wings  requires  special  mesh  moving  techniques 
as  well  as  techniques  to  represent  the  complex  deformation  patterns  of  the  wings.  In 
this  chapter  we  present  several  special  techniques  developed  by  the  Team  for  Advanced 
Flow  Simulation  and  Modeling  (T^AFSM)  to  address  these  and  other  computational 
challenges  related  to  modeling  of  a  flapping-wing  MAV.  For  completeness,  we  provide 
from  [1]  the  material  in  this  chapter. 

3.1  Time  representation 

Let  us  represent  time  t  G  (0,T)  with  order  NURBS  basis  functions,  (/9  = 
0, . . .  ,nct  —  1).  The  basis  functions  are  dehned  on  the  parametric  space  described 
by  the  open  knot  vector  {di, . . .  ,'dnkt},  where  net  and  nkt  are  the  number  of  control 
points  and  knots.  Then,  time  t  can  be  written  as 

t  =  (3.1) 

/3=0 

where  represents  the  temporal-control  point.  In  the  case  of  the  ST  formulation 
there  is  a  mesh  corresponding  to  each  temporal-control  point  tf.  As  originally  pro- 
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posed  in  [3],  we  use  a  strictly-increasing  mapping  function  0^(6*)  to  relate  the  element 
coordinate  for  a  knot  span  and  the  NURBS  parametric  space: 

^  ^  (1  —  0t(^))l^e+p+l  +  (1  +  Qt(^))'^e+p+2  ^2  2^ 

where  e  represents  the  element  index  (e  =  0, . . .  ,neit  —  1,  where  Ueit  is  the  number  of 
elements).  We  have  Uct  =  'n.kt  —  P  —  and  assuming  no  knot  multiplicity  inside  the 
knot  vector,  rieit  =  Ukt  —  2p  —  1.  The  element  shape  functions  are  dehned  as 

T^{Qt{e))  =  (i9) .  (3.3) 


In  the  time  interval  of  element  e,  we  represent  t  with  the  local  shape  functions: 

p+i 

(3.4) 

a=l 


3.1.1  Time  marching  problem 

Consider  a  set  of  NURBS  basis  functions  that  is  being  used  in  representing  the  data 
that  we  will  work  with.  Figure  3.1  shows  an  example.  Starting  from  this,  as  proposed 


CC 

4— < 

CC 

Q 


0.0  1.0  2.0  3.0  4.0 

Time  t 


Figure  3.1:  Data  represented  with  NURBS.  The  data  and  control  variables  (top). 
The  basis  functions  corresponding  to  each  control  variables  (bottom). 
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in  [3,  10],  we  can  form  a  new  basis  set  by  knot  insertion  with  the  objective  that  all 
elements  become  patches  as  shown  in  Figure  3.2.  After  that,  we  can  use  this  basis 
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Figure  3.2:  The  data  represented  with  basis  functions  after  the  knot  insertion.  The 
data  and  control  variables  (top).  The  basis  functions  corresponding  to  each  control 
variables  (bottom). 


set  in  our  ST  computation,  where  the  knot  spans  would  be  the  time  intervals  of  the 
ST  slabs,  and  represent  the  data  exactly.  An  alternative  process,  which  was  also 
proposed  in  [3,  10],  has  the  same  functionality,  but  without  the  need  to  explicitly 
represent  the  data  we  are  working  with  using  the  new  basis  set.  In  that  process, 
we  separately  form  a  new  basis  set  where  each  element  is  a  patch,  and  the  data  is 
represented  in  the  formulation  in  terms  of  its  own  basis  set.  In  general,  the  basis 
set  that  we  separately  form  does  not  need  to  be  the  same  as  the  one  that  could  be 
obtained  by  the  process  described  earlier,  but  if  it  is,  then  the  two  processes  result 
in  equivalent  solution  methods.  Figure  3.3  shows  the  new  basis  functions  that  we 
separately  form. 

Since  we  need  to  work  with  different  basis  sets,  we  map  one  parametric  space  to 
the  other  through  physical  time  t.  With  the  function  dehned  by  Eq.  (3.1),  time  t  can 
be  obtained  from  the  parametric  space 'd.  Here  we  consider  the  inverse  mapping;  i.e., 
t  ^ 'd.  Finding  the  parametric  space  coordinate  was  proposed  in  [3,  10]  as  follows: 
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Figure  3.3:  The  data  and  control  variables  (top).  The  basis  functions  that  we  sepa¬ 
rately  form  for  a  given  interval  for  the  ST  computation  (bottom).  To  integrate  over 
the  interval,  in  the  NURBS  representation  of  the  data  we  need  to  search  for  the  cor¬ 
responding  element  and  parametric  coordinate  for  the  time  t(d®)  of  each  quadrature 
point and  interpolate  the  value  from  the  data. 


1.  Find  the  element  e  that  is  represented  by  the  knot  span  ('dg+p+i, 'de+p+2)-  The 
process  requires  only  time  values  at  each  element  boundary,  and  we  can  quickly 
obtain  the  element  index  e  by  using  a  binary  search  technique. 


2.  Calculate  6  for  a  given  t  by  using  Newton-Raphson  iterations  as  follows: 


=  »'-(*- ye.  («■)))  I  ^ 


-1 


(3.5) 


where  “i”  is  the  iteration  counter,  t  (0*  (0*))  can  be  calculated  from  Eq.  (3.4), 
and 


dt 

de 


p+i 

E 

Q  =  1 


dT“ 

.e+a-l 


d0* 


d0* 


et(ep 


d9 


6» 


We  use  as  the  initial  guess  9^  =  0. 


(3.6) 


3.  Compute  -9  from  Eq.  (3.2). 
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3.1.2  Design  of  temporal  NURBS  basis  functions 

In  the  previous  section  we  described  how  to  hnd  the  parametric  space  value  cor¬ 
responding  to  physical  time.  Here  we  describe  from  [3,  10]  some  specihc  temporal 
representations. 

For  implementation  convenience  and  computational  efficiency,  we  restrict  the  time 
interval  of  the  ST  slab  such  that  the  time  interval  does  not  step  over  a  time  corre¬ 
sponding  to  a  temporal  knot  in  the  basis  set  used  for  representing  the  data  or  mesh. 
Thus  the  supporting  set  of  meshes  for  each  ST  slabs  consists  of  only  specihc  p  -|-  1 
meshes,  where  p  is  the  order  of  basis  used  for  representing  the  data  or  mesh.  Because 
of  that  requirement,  a  uniform  element  size,  i.e.  t('de+p+2)  —  t(de+p+i)  =  At,  where 
At  =  — ,  is  convenient.  Moreover,  we  might  have  the  following  requirement: 


dt  At 
d0  ^  T‘ 


(3.7) 


In  the  case  of  B-spline  basis  functions,  for  the  identity  mapping  0t(6')  =  6,  we  can 
satisfy  the  condition  expressed  by  Eq.  (3.7)  by  selecting  the  control  points  as  follows: 


+ 


•d/j+p+i  —  'd/3+i„ 

P  -  ^l) 


for  /d  =  1, . . . ,  Uct  —  1  and  =  0. 


(3.8) 


3.1.3  Approximation  in  time 

Let  Xa  sampling  values  of  a  time- varying  spatial  position  vector  Xyi  at  sampling 

times  (s  =  0,  •  ■  ■  ,  Ugp  —  1,  where  Ugp  is  the  number  of  sampling  points).  For 
example,  Xyi  could  be  the  position  vector  for  spatial  node  A,  or  it  could  be  the 
position  vector  for  a  point  on  a  surface  geometry  extracted  from  video  data.  For 
each  A,  as  proposed  in  [3,  10],  we  represent  the  path  corresponding  to  the  sampling 
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points  with  NURBS.  This  serves  two  purposes:  bringing  smoothness  to  the  temporal 
representation,  and  better  representation  accuracy  for  less  control  points.  First,  we 
form  a  linear  hnite  element  mesh  in  time,  consisting  of  two-node  elements.  Then,  we 
use  least-squares  projection  to  translate  that  to  NURBS  representation: 

[  R'X  ■  (x^  -  Xa]  =  0,  (3.9) 

Jo 

where  and  are  the  test  function  and  NURBS  representation  in  time,  and  x\  is 
the  linear  representation  in  time.  Thus,  we  obtain  the  control  point  corresponding 
to  each  time  control  point  t^.  Figure  3.4  is  an  example. 


Figure  3.4:  NURBS  representation  for  a  time-varying  spatial  position  vector.  The 
circles  are  the  spatial  position  vector  at  each  sampling  time.  The  squares  are  the 
temporal-control  points  and  the  smooth  curve  is  represented  by  them. 

Remark  5  This  is  a  simple  projection.  However,  the  concept  is  applicable  to  more 
complicated  formulations  to  obtain  smoother  motion. 

3.2  Simple-Shape  Deformation  Model  (SSDM) 

Suppose  we  want  to  track  the  motion/deformation  of  an  object  with  surface  shape  that 
is  too  complex  to  track  in  full  detail,  giving  us  just  the  option  of  tracking  only  a  hnite 
number  of  points  belonging  to  this  complex  shape.  In  the  simple-shape  deformation 
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Figure  3.5:  SSDM.  Complex  shape  is  shaded.  Circles  are  tracked  points.  SS  is 
represented  by  squares  (control  points). 

model  (SSDM),  we  assume  that  those  tracked  points  are  associated  with  a  simple 
shape  (SS)  instead  of  the  actual,  complex  shape.  NURBS  is  used  for  the  spatial 
representation  of  the  SS.  We  note  that  the  SS  is  larger  than  the  complex  shape.  The 
complex  shape  can  be  represented  by  hnite  elements  or  NURBS. 

Starting  with  the  reference  conhguration,  the  SS,  the  complex  shape  and  the 
tracked  points  all  can  be  seen  in  a  common  parametric  space  (Figure  3.5).  The 
complex  shape  can  be  represented  by  hnite  elements  or  NURBS.  Control  points  of 
the  SS  at  different  times  during  tracking  are  determined  by  a  least-squares  ht.  The  ht 
minimizes  the  difference  between  the  positions  on  the  SS  (with  respect  to  the  reference 
conhguration)  of  the  tracked  points  and  the  positions  of  the  actual  tracked  points.  The 
complex  shape  at  a  given  temporal-control  point  is  determined  by  interpolation  from 
the  parametric  space  in  the  case  of  the  hnite  element  representation,  and  by  least- 
square  projection  in  the  case  of  NURBS  representation.  The  least-squares  integration 
is  over  the  parametric  space  of  the  complex  shape,  and  we  minimize,  with  respect  to 
the  control  points  of  the  complex  shape,  the  diherence  between  the  complex-shape 
and  simple-shape  representations. 
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In  the  full  ST  representation,  the  method  we  described  above  is  applied  to  temporal- 
control  values  that  are  determined  as  described  in  Section  3.1.3,  instead  of  the  actual 
physical  locations. 

3.3  Mesh  update  techniques  in  the  temporal  NURBS 
representation 

3.3.1  Mesh  computation  and  representation 

Given  the  surface  mesh,  we  compute  the  volume  mesh  using  the  mesh  moving  tech¬ 
nique  we  have  been  using  [52].  Here,  as  proposed  in  [3],  we  apply  this  technique 
to  computing  the  meshes  that  will  serve  as  temporal-control  points.  This  allows 
us  to  do  mesh  computations  with  longer  time  in  between,  but  get  the  mesh-related 
information,  such  as  the  coordinates  and  their  time  derivatives,  from  the  temporal 
representation  whenever  we  need.  Obviously  this  also  reduces  the  storage  amount 
and  access  associated  with  the  meshes.  However,  because  of  the  longer  time  be¬ 
tween  the  control  meshes,  linear  interpolation  of  the  surfaces  between  control  points 
in  time  might  be  needed  in  computing  those  meshes  with  the  mesh  moving  technique 
mentioned. 

Remark  6  We  note  that  getting  the  meshes  used  in  the  computations  from  the  tem¬ 
poral  representation  can  be  done  independent  of  which  time  direction  was  used  in  com¬ 
puting  the  control  meshes.  In  other  words,  it  does  not  matter  if  the  control  meshes 
were  computed  while  the  wings  were  flapping  forward  or  backward  in  time. 

3.3.2  Remeshing  techniques 

In  many  computations  remeshing  becomes  unavoidable.  Two  choices  were  proposed 
in  [3].  To  explain  those  two  choices,  let  us  assume  that  when  we  try  to  move  from 
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control  mesh  to  we  find  the  quality  of  to  be  less  than  desirable. 

In  the  first  choice,  which  is  called  “trimming,”  we  remesh  going  back  to  . 

Then  whenever  our  solution  process  needs  a  mesh,  depending  on  the  time,  we  use  the 
control  meshes  belonging  to  either  only  the  un-remeshed  set  or  only  the  remeshed  set 
(Figure  3.6). 

In  the  second  choice,  we  perform  knot  insertion  p  times  in  the  temporal  represen¬ 
tation  of  the  surface  at  the  right-most  knot  before  the  maximum  value  of  the  basis 
function  corresponding  to  making  that  knot  a  new  patch  boundary.  Then  we 

do  the  mesh  moving  computation  for  the  control  meshes  associated  with  the  newly- 
defined  basis  functions,  not  only  the  one  at  the  new  patch  boundary,  but  also  going 
back  (p  —  1)  basis  functions  (Figure  3.7). 

We  use  the  second  choice  in  computations,  because  we  believe  that  in  many  cases 
the  need  for  remeshing  is  generated  by  a  topological  change,  which  we  can  avoid  going 
over  with  a  large  step  if  we  use  the  knot  insertion  process. 

3.4  Fluid  mechanics  computation  with  temporal 
NURBS  mesh 

We  solve  the  fluid  dynamics  equations  with  the  DSD/SST  formulation.  Here  we 
describe  two  techniques  related  to  the  moving-mesh  problem. 

3.4.1  No-slip  condition  on  a  prescribed  boundary 

Suppose  we  have  a  prescribed  mesh  motion,  and  no-slip  conditions  on  part  of  the 
boundary  of  that  mesh.  Those  Dirichlet  conditions  can  be  obtained  from  the  mesh 
boundary  motion. 

Prior  to  solving  the  equations  using  a  ST  slab  Qn,  we  use  a  least-squares  projection 
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Figure  3.6:  Remeshing  and  trimming  NURBS.  A  set  of  un-remeshed  meshes  (top). 
A  set  of  remeshed  meshes  (middle).  Common  basis  functions  (bottom). 
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Figure  3.7:  Remeshing  with  knot  insertion.  For  the  set  of  un-remeshed  meshes,  there 
are  p  newly-dehned  basis  functions  and  the  corresponding  control  points  are  marked 
“New.”  We  carry  out  the  mesh  moving  computations  for  those  meshes. 
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for  each  prescribed  node  A  as  follows: 

/tn  +  l  /  \ 

^  j  dt  =  0,  (3.10) 

where  represents  the  test  function,  is  represented  by  temporal-control  veloci¬ 
ties  (unknown)  and  the  corresponding  basis  functions  in  time,  and  the  mesh  velocity 
is  obtained  by  the  derivative  of  the  mesh  displacement,  which  is  also  represented  by 
temporal-control  positions  and  their  basis  functions.  We  note  that  at  time 
approaching  from  below  and  above  might  be  different. 


3.4.2  Starting  condition 

Starting  a  fluid  dynamics  computation  is  not  always  easy,  especially  in  the  presence 
of  moving  boundaries.  For  example,  pre-FSI  computation  techniques  were  developed 
in  [113,  89,  86,  72]  to  build  a  starting  condition  for  FSI  computations.  Here  we 
describe  the  techniques  we  use  as  precomputation  sequence  for  flow  computations 
with  prescribed  boundary  motion. 

Suppose  we  want  to  compute  with  a  mesh  temporally  represented  with  NURBS 
(M°,  M^,  ■  ■  ■).  We  define  the  temporal-control  value  to  be  used  interchangeably 
with  the  temporal-control  mesh  M^. 


With  quadratic  NURBS  representation 

We  describe  this  option  from  [3].  Figure  3.8  is  an  example  of  mesh  representation 
with  quadratic  NURBS  in  time.  We  generate  two  additional  meshes  as  follows: 


M-i  = 


-  aM} 


a 


=  M, 


-1 


(3.11) 


(3.12) 
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Figure  3.8:  Mesh  representation  for  the  starting  condition  with  quadratic  NURBS  in 
time. 


where  0  <  a  <  1  is  an  extrapolation  parameter.  Mesh  ^  is  an  extrapolation.  The 
corresponding  temporal-control  point  for  M~^  is 


t 


-1 

C 


tp  - 
1  —  a  ’ 


(3.13) 


with  the  only  requirement  being  t~^  =  tg  <  t~^,  which  determines  the  length  of  the 
precomputation.  For  the  computations  reported  in  this  paper,  in  temporal  represen¬ 
tation  of  the  mesh,  as  NURBS  basis  functions,  we  use  quadratic  B-spline  functions 
dehned  by  the  knot  vector  {0,  0,  0, 1, 1, 1}.  Based  on  the  preliminary  computations, 
using  even  higher-order  NURBS  basis  functions  was  proposed  in  [3]  so  that  the  ac¬ 
celeration  is  continuous.  We  explain  this  in  the  next  subsection. 


With  cubic  NURBS  representation 

Figure  3.9  is  an  example  of  mesh  representation  with  cubic  NURBS  in  time.  Here 
we  construct  a  starting- condition  temporal  patch  with  U^-continuity  with  the  next 


53 


35 


M, 


-3 


^0. 


M: 


H- 


Mi 


^3  <  t  <  0  0  <  t 

Figure  3.9:  Mesh  representation  for  the  starting  condition  with  cubic  NURBS  in  time. 


patch: 


lim  X 

t-5-O- 


lim 

t — >^0 


dx 

dt 


.U2 

t^o-  dt^ 


lim  X, 

t^o+ 


lim 

t-s>o+  dt 

d^x 

lim 

t-s-o+  dt"* 


(3.14) 

(3.15) 

(3.16) 


where  x  represents  the  mesh  position  dehned  with  temporal-control  meshes,  and  we 
want  the  mesh  velocity  and  acceleration  to  be  zero  at  the  beginning  of  that  temporal 
patch: 


dx 

dt 

d^x 


dt^ 


t  =  ts 


0, 

0. 


(3.17) 

(3.18) 


Remark  7  Because  this  condition  provides  zero  velocity  and  acceleration,  we  can 
compute  with  a  static  mesh  prior  to  the  computation  in  this  temporal  patch  to  develop 
the  flow  field. 
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To  achieve  above  conditions  we  use  cubic  B-spline  functions  defined  by  the  knot 
vector  {0,  0,  0,  0,  1, 1, 1, 1}.  Because  the  last  temporal-control  mesh  for  the  starting- 

condition  patch  is  the  condition  given  by  Eq.  (3.14)  is  automatically  satished. 
The  conditions  given  by  Eqs.  (3.17)  and  (3.18)  can  be  satished  by  setting  = 
M~^  =  M~‘^,  with  any  set  of  corresponding  control  points  for  time:  t~‘^  =  t^<  < 

t~‘^.  However,  satisfying  the  conditions  given  by  Eqs.  (3.15)  and  (3.16)  is  not  trivial. 

The  following  is  an  example  of  using  cubic  B-splines.  The  derivatives  of  the  mesh 
position  at  t  — )■  0^  are  as  follows: 


dx 


dd 


d^x 


dd2 


(3.19) 

(3.20) 


where  d  is  parametric  coordinate^  dehned  on  the  hrst  patch.  Similarly,  the  derivatives 
of  the  time  are  as  follows: 


dt 

dd 

d^t 


dd2 


_ 7 _ 

^5-^2 

6  (tl-tl  _  tl  \ 

1^5  -1^3  We  -^3  -^^2/  ’ 


(3.21) 

(3.22) 


Thus,  the  mesh  velocity  and  acceleration  can  be  determined  as  follows: 


dx 
hm  — 

t-5>o+  dt 


,  d^x 
hm  — - 
t^o+  dt^ 


dtyVd^x  dxd^tx 
ddy 


(3.23) 

(3.24) 


with  Eqs  (3.20)-(3.22).  We  note  that  if  the  temporal-control  points  are  given  by 
Eq.  (3.8),  the  second  term  of  Eq.  (3.24)  is  zero. 


^We  note  that  the  index  of  the  knot  vector  starts  from  1. 
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Similarly,  the  derivatives  at  t  — )■  0  are  as  follows: 


dx 

lim  —  =  — 

t-5>o-  at 

d^x  1 
lim  ; 

t^o-  dt^  d{t-^y 


-  3x,  +  2x,  + 


tc  -  3t, 
+-1 


-1 


(3.25) 


(3.26) 


With  that,  there  are  many  ways  of  satisfying  the  conditions  given  by  Eqs.  (3.15) 
and  (3.16).  Here  we  set  t~^  =  —Atg.  We  also  dehne  the  control  points  for  the  time  in 
snch  a  way  that  the  part  of  Eq.  (3.26)  related  to  ^  at  t  — )■  0“  is  zero;  i.e., 

Thus,  we  can  set  x~^  and  x“^  as 


Atg 


3 


?)  -  Zli?  lim 

^  dt2  J 


(3.27) 

(3.28) 


Furthermore  we  set  t~^  =  —5Atg  and  tg  =  —6Atg  to  have  an  overall  linear  relationship 
between  the  time  and  the  parametric  space.  We  leave  Atg  as  a  free  parameter  that 
we  can  adjust  for  various  purposes,  such  as  having  smaller  acceleration  and  smaller 
displacement  from  M°. 

Remark  8  In  addition  to  using  NURBS  basis  functions  for  the  temporal  representa¬ 
tion  of  the  motion  and  deformation  of  the  solid  surfaces  and  volume  meshes  computed, 
we  can  use  NURBS  as  temporal  basis  functions  in  our  ST  computations.  With  that, 
we  would  have  temporal  advantages  similar  to  the  spatial  advantages  one  would  have 
by  using  NURBS  as  spatial  basis  functions. 
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Chapter  4 


Flapping-motion  and  geometry 
representation 


In  this  chapter,  we  describe  from  [1]  the  techniques  used  for  modeling  the  body 
and  wings  of  the  MAV  and  the  techniques  used  for  representing  the  motion  and 
deformation  of  the  flapping  wings,  with  the  data  coming  from  the  actual  locust. 

4.1  Wing  and  body  shape 

Based  on  a  digital  image  of  a  pair  of  locust  forewing  (FW)  and  hindwing  (HW),  we 
construct  the  surface  geometries  of  the  FW  and  HW  using  a  NURBS-based  design 
software  (see  Figure  4.1).  The  FW  and  HW  are  modeled  as  a  single  NURBS  patch 
with  degeneration  at  the  wingtip.  There  are  21  control  points  for  the  FW,  and  51 
for  the  HW  (see  Figure  4.2).  in  [1].  The  HWs  for  the  MAV  are  attached  to  the  body 
through  the  entire  chord  length  while  the  HWs  for  the  locust  computation  presented 
in  [1]  are  not  (see  Figure  4.3). 

The  MAV  body  is  a  simplihed  and  streamlined  version  of  the  locust  body.  The 
body  is  symmetric  along  the  sagittal  plane.  The  body  is  made  up  of  20  NURBS 
patches,  10  patches  in  each  sagittal  half  (this  is  the  minimum  number  needed  to 
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Figure  4.1:  Wing  model  construction. 


Figure  4.2:  MAV  FW  and  HW  surfaces  represented  by  NURBS  and  the  corresponding 
control  mesh. 
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Figure  4.3:  Locust  body  and  wings,  see  [1],  (left)  and  MAV  body  and  wings  (right), 
attach  the  wings  to  the  body).  We  shape  the  patches  as  shown  in  Figure  4.4.  The 


Figure  4.4:  MAV  body  patches. 

wing  spatial-control  meshes  are  deformed  to  the  various  positions  in  the  flapping 
motion  as  we  will  describe  in  Section  4.2. 


4.2  Flapping- mot  ion  representation 

The  wing  motions  we  prescribe  need  to  accurately  represent  the  wing  motion  and 
deformation  patterns  during  flight.  Here  we  face  the  challenge  of  reconciling  data 
acquired  through  photogrammetry  with  data  that  is  suitable  as  an  input  for  compu- 
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tational  analysis.  In  this  compntation,  we  nse  a  data  set  that  is  more  recent  than  the 
one  reported  in  [3].  Both  data  sets  were  provided  by  onr  collaborators  at  BCM.  The 
location  of  the  tracking  points  provided  in  this  data  set  can  be  seen  in  Fignre  4.5. 
By  averaging  the  valnes  corresponding  to  each  other  across  the  symmetry  plane,  we 
create  a  new  data  set  of  tracking  points  for  one  side.  With  the  help  of  these  tracking 
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Fignre  4.5:  Tracking  points  in  data  set  provided  by  BCM. 

points,  we  generate  additional  “tracking  points”  to  more  closely  represent  locnst  flight 
characteristics  observed  in  the  video.  After  the  body  motion  is  remove,  the  motion 
of  the  wings  is  reflected  across  the  sagittal  plane  of  the  MAV  to  create  symmetric 
flapping  motion.  The  total  nnmber  of  tracking  points  nsed  in  the  compntations  is  68. 
We  then  use  this  symmetric  data  as  the  input  for  the  least-squares  projection  used 
for  the  temporal  representation.  Next,  we  temporally  interpolate  our  representative 
data  set.  For  each  tracking  point,  we  apply  a  temporal  NURBS  representation  as 
discussed  in  Section  3.1.2  and  3.1.3.  Higher-order  basis  functions  should  be  used  in 
temporal  interpolation  to  achieve  a  continuous  acceleration.  We  use  cubic  B-spline 
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functions  for  temporal  interpolation.  Motion  with  cubic  basis  functions  maintains  C^- 
continuity  across  knots.  With  this,  acceleration  will  be  continuous  across  temporal 
knots. 

4.2.1  Generating  a  periodic  data  set 

While  flapping  motion  kinematics  is  inherently  cyclical,  it  is  not  necessarily  periodic. 
For  example,  the  vertical  position  at  the  start  of  the  first  downstroke  does  not  corre¬ 
spond  to  exactly  the  same  position  for  subsequent  strokes.  This  adds  to  the  difficulty 
of  computing  flapping-flight  aerodynamics.  It  is  beneficial  to  have  a  periodic  data  set 
for  a  single  flap  cycle,  where  the  first  and  last  points  of  the  cycle  are  colocated  and 
have  the  desired  continuity.  Thus,  a  single  set  of  deformed  meshes  can  be  appended 
to  produce  as  many  flapping  cycles  as  are  reqnired.  Because  of  the  periodicity  of  the 
motion,  we  can  compute  a  desire  number  of  flapping  cycles  that  we  think  would  be 
sufficient  to  reach  a  solution  that  we  would  consider  periodic. 

To  obtain  a  periodic  data  set,  after  the  least-squares  projection,  we  extract  one 
flapping  cycle  by  averaging  and  inserting  knots  at  the  top  of  the  HW  cycle  (end  of 
the  npstroke  and  beginning  of  the  downstroke).  To  maintain  continnity,  the  control 
points  corresponding  to  the  knot  at  the  beginning  of  the  flap  cycle  are  colocated  with 
the  control  points  corresponding  to  the  knot  at  the  end  of  the  flap  cycle  (for  a  cubic 
B-spline,  three  control  points  correspond  to  a  given  knot).  To  obtain  snch  repetition 
we  average  those  control  points.  Finally,  we  insert  knots  to  extract  a  single  cycle.  We 
show  the  averaging  process  for  a  cnbic  B-spline  in  Figure  4.6. 

4.2.2  Motion  of  the  wings 

We  spatially  represent  the  tracking  points  at  each  temporal-control  point.  Spatial 
interpolation  is  accomplished  using  the  SSDM  described  in  Section  3.2.  The  FW  and 
HW  SS  consist  of  6  and  9  control  points,  respectively. 
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(a)  Flap  cycle  of  interest  before  averaging  con¬ 
trol  points  (left)  and  enlarged  view  of  top  of  cycle 
(right).  The  control  points  we  average  are  high¬ 
lighted. 


(b)  Flap  cycle  of  interest  after  averaging  control 
points  (left)  and  enlarged  view  of  top  of  cycle 
(right). 


Figure  4.6:  Process  used  for  building  a  periodic  flapping  cycle. 
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Chapter  5 


Base  MAV  computation 


In  this  chapter,  we  describe  the  MAV  computations  where  the  wing  deformation 
patterns  come  from  the  actual  locust.  These  deformation  patterns  and  the  motion 
of  the  wings  are  prescribed  as  described  in  Chapter  4.  In  the  rest  of  the  chapter, 
from  [82] ,  we  describe  the  setup  for  the  computation,  mainly  the  surface  and  volume 
meshes  and  the  mesh  update  technique,  and  report  the  aerodynamic  forces  generated. 

5.1  Surface  and  volume  meshes 

After  the  wing  spatial-control  surfaces  are  deformed  to  the  various  positions  in  the 
flapping  motion  as  described  in  Section  4.2.2,  we  discretize  the  surfaces  at  each 
temporal-control  point  with  triangular  elements.  The  triangular  surface  mesh  used 
in  the  computation  is  shown  in  Figure  5.1.  We  note  that  both  FW  and  HW  have  a 
hnite  thickness  of  1%  of  the  FW  root  chord,  which  tapers  to  zero  thickness  at  the 
wing  edges. 

We  generate  the  tetrahedral  element  volume  mesh  with  the  following  steps.  First, 
we  generate  a  one-layer  rehnement  region  near  the  wing  surface.  The  element  height 
of  this  layer  is  10%  of  the  FW  root  chord.  In  addition,  we  have  a  cylindrical  region 
of  increased  element  rehnement  around  the  MAV.  We  also  specify  increased  volume- 
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Figure  5.1:  MAV  wing  and  body  surface  meshes  with  triangular  elements. 

mesh  rehnement  in  the  region  between  the  FW  and  HW.  This  rehnement  reduces 
the  mesh  distortion  due  to  the  staggered  motion  of  FW  and  HW.  The  volume  mesh 
within  this  cylindrical  region  is  then  generated  using  an  automatic  mesh  generator. 
We  rotate  the  cylindrical  region  to  an  angle  representing  the  approximate  free-flight 
body  angle  of  the  locust  and  compute  mesh  motion  only  in  this  cylindrical  region 
as  discussed  in  Section  4.2.  A  volume  mesh  between  the  cylindrical  region  and  the 
external  domain  boundaries  is  then  generated  with  an  automatic  mesh  generator. 
This  region  is  hxed  for  the  entire  computation;  therefore,  it  is  used  for  all  temporal 
patches.  The  volume  mesh  boundaries  and  cylindrical  rehnement  region  are  shown 
in  Figure  5.2. 

5.2  Mesh  update  technique 

Because  of  the  topological  changes  caused  by  the  two-wing  motion,  we  cannot  use 
the  same  connectivity  during  the  entire  computation.  We  know  the  wing  motion  in 
advance,  so  we  use  a  “prepared”  remeshing  technique. 
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Figure  5.2:  Surface  meshes  on  some  of  the  external  computational  boundaries  and 
on  the  outer  surface  of  the  inner,  cylindrical  mesh  region,  which  has  been  rotated  to 
account  for  in-flight  body  angle. 


65 


47 


We  have  a  surface  modeling  using  a  temporal  NURBS  representation.  First  we 
decide  the  remeshing  instants.  Then,  we  use  knot  insertion  to  split  the  flapping  cycle 
into  patches  as  seen  in  Section  3.3.2.  Those  patches  are  called  “temporal  patches.” 
The  number  of  nodes  and  elements  in  the  total  volume  mesh  varies  between  each 
temporal  patch.  The  average  number  of  nodes  and  elements  are  0.35  million  and  2.1 
million  (see  Table  5.1).  Due  to  the  relatively  large  change  in  deformation  between 
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Temporal 
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Patch 

Point 

of  Nodes 

Elements 

1 

4 

347,312 

2,072,463 

2 

3 

342,501 

2,043,814 

3 

3 

334,840 

1,998,235 

4 

3 

385,900 

2,303,385 

Table  5.1;  Summary  of  computational  setup.  In  each  temporal  patch  we  indicate  the 
number  of  temporal  control  points  and  at  which  point  we  generate  the  tetrahedral 
volume  mesh  with  the  number  of  nodes  and  elements  shown. 


each  temporal-control  point,  we  use  sub  iterations  for  the  mesh  computation  to  divide 
the  steps  between  temporal-control  points  into  20  smaller  steps.  We  move  the  mesh, 
which  corresponds  to  the  middle  control  point,  backward  and  forward  through  each 
smaller  step  using  1,500  GMRES  [114]  iterations  (this  is  the  technique  described  in 
Section  3.3.1). 


5.3  Computation  of  the  base  case 

5.3.1  Computational  conditions 


Figure  5.3  shows  the  length  scales  involved  in  the  model  used  in  the  computations. 
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Figure  5.3:  Length  scales  involved  in  the  model  used  in  the  computations. 

The  air  density  and  kinematic  viscosity  are  1.2251  kg/m^  and  1.4606x10“^  m^/s. 
The  period  of  flapping  is  T  =  0.047  s.  Prior  to  the  beginning  of  the  prescribed 
flapping  motion,  we  compute  400  time  steps  to  develop  the  flow  held.  Over  the 
hrst  300  time  steps  of  this  computation,  we  use  a  Cosine  form  to  smoothly  increase 
the  inhow  velocity  from  0.0  to  2.4  m/s,  which  represents  the  average  wind  tunnel 
velocity  used  for  the  empirical  data  collection.  For  the  last  100  time  steps  we  keep  the 
velocity  steady  at  2.4  m/s.  In  this  how-development  computation,  the  time  step  size 
is  1.71x10“^  s,  with  4  nonlinear  iterations  per  time  step.  We  use  the  ST-SUPS  and 
ST-VMS  techniques  for  the  hrst  two  and  last  two  nonlinear  iterations,  respectively, 
with  the  stabilization  parameters  as  given  by  Eqs.  (7)-(ll)  in  [8]  for  Tm  =  Tsupg  (which 
are  also  in  Section  2.3)  and  Eqs.  (35)-(37)  in  [3]  for  Vc-  The  stabilization  parameters 
are  calculated  after  the  predictor  step  and  after  the  hrst  two  nonlinear  iterations.  The 
number  of  GMRES  iterations  for  the  nonlinear  iterations  are  30,  60,  120,  and  180. 

For  the  computation  with  happing,  we  use  25  ST  slabs  (with  linear  basis  functions) 
for  each  of  the  knot  spans  in  the  temporal  representation  of  the  mesh,  and  that  gives 
us  the  time-step  size  of  1.71x10“^  s.  The  rest  of  the  computational  parameters  are 
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the  same  as  those  above. 

Remark  9  We  have  seldom  observed  close-to-zero  or  negative  diagonal  terms  when 
the  time-step  size  is  large.  This  occurs  when  the  fine-scale  velocity  is  much  larger 
than  the  coarse-scale  (discrete)  velocity.  We  believe  this  is  due  to  the  predictor,  which 
assumes  all  velocities  are  the  same  from  the  previous  time  step  except  those  with 
Dirichlet  conditions. 

We  partition  this  mesh  to  enhance  parallel  efficiency.  Mesh  partitioning  is  based 
on  the  METIS  [115]  algorithm.  We  partition  it  into  128  parts  and  nse  128  parallel 
processes. 


5.4  Results 


We  compute  the  problem  for  three  complete  flapping  cycles.  We  display  results  for  the 
third  cycle.  The  lift  and  thrust  (pressure  component  only)  are  shown  in  Figures  5.4- 
5.6  for  the  MAV  base  case.  Figures  5. 7-5. 8  show  the  pressure  on  the  wing  surfaces 


Figure  5.4:  Total  lift  and  thrust  generated  over  one  cycle.  Positive  value  indicates 
that  thrust  exceeds  drag. 
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Figure  5.6:  Lift  and  thrust  generated  on  the  right  HW. 
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(relative  to  the  free-stream  pressure)  at  different  instants  during  the  third  flapping 
cycle.  Figure  5.9  shows  the  vorticity  magnitude  at  different  instants. 


Figure  5.7:  Surface  pressure  in  Pa  (relative  to  the  free-stream  pressure)  at  the  hrst 
four  of  the  eight  equally-spaced  points  during  the  third  flapping  cycle  (top  view  on 
left,  bottom  view  on  right). 
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Figure  5.8:  Surface  pressure  in  Pa  (relative  to  the  free-stream  pressure)  at  the  last 
four  of  the  eight  equally-spaced  points  during  the  third  flapping  cycle  (top  view  on 
left,  bottom  view  on  right). 
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Figure  5.9:  Volume  rendering  of  the  vorticity  magnitude  for  the  eight  equally-spaced 
points  during  the  third  flapping  cycle  (left  to  right  and  then  top  to  bottom). 


73 


Chapter  6 


SCFSI  computation 


In  this  chapter,  we  describe  from  [2]  the  SCFSI  computations  for  the  MAV.  We  do  the 
structural  mechanics  computation  only  for  the  HW,  as  it  generates  most  of  the  lift 
and  thrust.  We  use  the  Kirchhoff-Love  shell  model.  The  FW  has  prescribed  motion 
and  deformation  as  described  in  Chapter  5. 

6.1  Setup  and  computations 

To  make  the  structural  mechanics  computations  less  challenging,  instead  of  the  com¬ 
plex  shape  for  the  HW,  we  use  the  SS.  We  use  a  higher  spatial  rehnement  for  the 
SS  than  the  one  in  [82].  It  has  16  control  pointed  and  is  represented  with  quadratic 
NURBS  in  both  directions  (there  were  9  control  points  in  [82]).  The  SS  for  the  FW  has 
the  same  number  of  control  points  (6)  as  in  [82].  Figure  6.1  shows  the  spatial-control 
mesh  for  the  HW  SS. 

We  start  with  a  structural  mechanics  computation  with  no  fluid  mechanics  forces 
applied.  We  prescribe  the  leading-edge  motion  of  the  wing.  We  can  use  either  a 
desired  leading-edge  motion  or  the  original  leading-edge  motion  from  the  locust  data, 
which  is  what  we  do  here.  We  note  that  the  SS  leading-edge  data  is  represented 
temporally  by  cubic  NURBS  and  can  easily  be  altered  while  maintaining  the  smooth 
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Figure  6.1:  Spatial-control  mesh  for  the  HW  SS. 


motion  and  continuous  acceleration.  The  structural  mechanics  material  properties 
are  based  on  measured  values  of  the  actual  locust  wing  material  properties,  with  the 
values  altered  to  take  into  account  the  homogeneous  modeling  of  the  wing  rather 
than  using  the  actual  wing  morphology  [116].  The  Young’s  modulus  is  4.3  GPa,  the 
density  is  560  kg/m^,  and  the  thickness  is  1.1x10“^  m. 

In  the  spatial  discretization  with  quadratic  NURBS  in  each  direction  and  16  con¬ 
trol  points  (see  Figure  6.1),  we  have  3  quadrature  points  in  each  direction,  and  the 
mass  matrix  is  consistent.  In  the  temporal  discretization,  we  use  the  generalized-a 
method  [112].  The  parameters  we  use  with  the  method,  in  the  notation  of  [11],  are 
Om  =  2,  Of  =  1,  7  =  1.5,  and  (3  =  1.  These  values  result  in  a  spectral  radius  of  zero 
while  maintaining  second-order  accuracy  in  time.  The  acceleration  predictor  is 

a°+i  =  a„,  (6.1) 


which,  because  of  the  Newmark- velocity  formula,  leads  to  the  following  displacement 
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and  velocity  predictors: 


Yn+i  =  Yn  +  u„Z\t  +  (6.2) 

u°+i  =  u„  +  a„Z\t.  (6.3) 

The  time-step  size,  1.71x10“^  s,  is  the  same  as  that  of  the  fluid  mechanics  computa¬ 
tions,  and  the  number  of  nonlinear  iterations  at  each  time  step  is  7.  The  number  of 
GMRES  iterations  for  each  nonlinear  iteration  is  sufficiently  high  to  make  the  solu¬ 
tion  process  equivalent  to  a  direct  one.  The  preconditioner  is  nodal-block-diagonal. 
We  compute  three  full  flapping  cycles.  We  use  the  motion  and  deformation  from  the 
third  cycle  in  the  fluid  mechanics  computation.  This  computed  data  replaces  the 
original  tracking-point  data  for  the  locust  HW,  and  the  rest  of  the  setup  process, 
mesh  generation  and  moving,  and  computational  parameters  are  almost  the  same 
as  those  described  in  Section  5.3.  The  number  of  GMRES  iterations  for  the  non¬ 
linear  iterations  are  a  little  different:  30,  60,  120,  and  200,  and  the  preconditioner 
is  nodal-block-diagonal,  instead  of  diagonal.  Also,  we  split  the  flapping  cycle  into 
Eve  temporal  patches  instead  of  four,  because  the  increased  complexity  of  the  wing 
deformation  patterns  requires  one  additional  remesh.  As  shown  in  [82],  increasing  the 
number  of  temporal  patches  beyond  four  does  not  result  in  a  significant  change  in  the 
aerodynamic  forces.  Figure  6.2  shows  the  wing  motion  and  deformation  used  in  the 
fluid  mechanics  computation  of  the  first  SGFSI  iteration.  We  do  the  fluid  mechanics 
computation  for  three  flapping  cycles. 

In  the  second  SGFSI  iteration,  the  structural  mechanics  computation  is  carried  out 
with  the  fluid  mechanics  forces  applied  from  the  fluid  mechanics  computation  for  the 
third  flapping  cycle  of  the  first  iteration.  We  calculate  the  pressure  difference  between 
the  top  and  bottom  surface  for  each  node  of  the  wing  at  each  time  step.  Figure  6.3 
shows  this  net  force  distribution  at  different  instants  during  the  third  flapping  cycle. 
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Figure  6.2:  MAV  with  the  HW  deformation  extracted  from  the  structural  mechanics 
computation,  seen  at  eight  equally-spaced  instants  during  the  third  flapping  cycle. 
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We  project  this  force  distribution  to  a  temporal  representation  with  cubic  NURBS. 
Then  we  project  that  from  the  complex  shape  used  in  the  fluid  mechanics  computation 
to  the  SS  used  in  the  structural  mechanics  computation.  We  can  discretize  that  in 
time  to  be  used  in  the  structural  mechanics  computations.  Temporal  representation 
of  the  force  with  NURBS  allows  us  to  discretize  in  time  with  any  time-step  size, 
without  loosing  the  smoothness  of  the  force  distribution  data.  We  perform  four 
SCFSI  iterations  where  this  process  is  repeated. 

6.2  Results 

We  present  the  results  for  each  SCFSI  iteration  compared  to  the  base  case.  Figure  6.4 
shows  the  total  lift  and  thrust.  Figures  6.5  and  6.6  show  the  lift  and  thrust  generated 
by  the  HW  and  FW. 

Even  though  the  FW  motion  is  not  different  from  what  it  was  in  the  base  case, 
because  of  the  interaction  with  the  HW,  the  force  the  FW  experiences  is  different. 
We  note  that  there  is  not  much  change  from  the  third  iteration  to  the  fourth,  and 
not  even  from  the  second  iteration  to  the  third,  and  therefore  we  will  stop  the  SCFSI 
iterations  at  three  in  the  computations  presented  in  the  rest  of  the  paper. 
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Figure  6.3:  Net  force  distribution  (in  Pa)  at  eight  equally-spaced  instants  during  the 
third  flapping  cycle. 
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Chapter  7 


Membrane  model 


It  is  interesting  to  investigate  what  effect  the  structural  model  has  on  the  performance 
of  the  flapping- wing  MAV.  In  this  chapter,  from  [2],  we  describe  the  computations 
with  the  membrane  model.  This  is  conceivable  if  the  MAV  wings  are  made  of  very  thin 
flexible  material  (on  the  order  of  micrometers)  with  negligible  bending  stiffness.  The 
leading  edge,  with  its  motion  prescribed,  behaves  as  a  beam-like  structural  member. 
In  an  actual  flapping-wing  MAV,  such  a  design  is  quite  plausible. 

7.1  Setup  and  computations 

With  the  membrane  model  we  use  the  same  procedure  and  computational  settings  as 
with  the  shell  model  in  Chapter  6.  We  also  use  the  same  material  properties.  After 
three  cycles  of  structural  mechanics  computation  for  the  first  SCFSI  iteration,  we 
proceed  as  before  with  three  cycles  of  fluid  mechanics  computation.  Figures  7.1  and 

7.2  show,  for  the  shell  and  membrane  models,  the  wing  motion  and  deformation  used 
in  the  fluid  mechanics  computation  of  the  first  SCFSI  iteration.  As  can  be  seen,  the 
membrane  model  results  in  larger  deformations.  We  perform  three  SCFSI  iterations, 
and  present  the  results  from  that. 
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Figure  7.1:  Top  view  of  the  MAV  with  the  HW  deformation  extracted  from  the 
structural  mechanics  computation,  seen  for  the  shell  model  at  eight  equally-spaced 
instants  during  the  third  flapping  cycle.) 
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Figure  7.2:  Top  view  of  the  MAV  with  the  HW  deformation  extracted  from  the 
structural  mechanics  computation,  seen  for  the  membrane  model  at  eight  equally- 
spaced  instants  during  the  third  flapping  cycle.) 
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7.2  Results 

Figures  7.3  and  7.4  show  the  total  and  right  HW  lift  and  thrust  with  the  shell  and 
membrane  models.  The  average  total  lift  with  the  membrane  model  is  1.70  mN 
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Figure  7.3:  Total  lift  (top)  and  thrust  (bottom)  with  the  shell  and  membrane  models. 


compared  to  1.43  mN  with  the  shell  model.  The  average  thrust  is  —2.04  mN  with  the 
membrane  model  and  —0.56  mN  with  the  shell  model.  This  shows  a  small  increase  in 
lift  with  a  signihcant  decrease  in  thrust.  However,  it  is  interesting  to  note  that  while 
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Figure  7.4:  Right  HW  lift  (top)  and  thrust  (bottom)  with  the  shell  and  membrane 
models. 
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the  downstroke  with  shell  model  produces  more  lift,  the  upstroke  with  the  membrane 
model  produces  quite  a  bit  of  positive  lift  instead  of  negative  lift.  This  suggests  that 
an  optimal  flapping-wing  MAV  design  should  perhaps  have  a  shell-like  behavior  on 
the  downstroke  and  membrane-like  behavior  on  the  upstroke.  Such  a  wing  would 
also  more  accurately  mirror  the  actual  locust  wing  deformation  in  flight.  Adding 
structural  members  that  provide  stiffness  to  a  membrane  wing  on  the  downstroke  is 
a  possible  solution. 


Chapter  8 


Prescribed  body  motion 


This  chapter  is  from  [2].  To  fully  simulate  an  MAV  under  various  flying  conditions, 
we  want  to  be  able  to  not  only  change  the  wing  deformation  patterns  and  statically 
change  the  angle  of  attack,  but  also  have  the  ability  to  apply  any  prescribed  body 
motion  to  the  MAV.  We  can  add  heave,  sway,  snrge,  roll,  pitch,  and  yaw,  and  in¬ 
vestigate  their  effect  on  the  forces  generated  by  the  MAV,  assnming  that  the  MAV 
somehow  performs  snch  manenvers.  Thanks  to  the  rehnement  cylinder,  we  can  add 
all  of  these  without  impacting  the  mesh  quality  around  the  MAV  wings.  All  the  mesh 
motion  is  taken  up  by  the  external  region  where  the  elements  can  nndergo  large  defor¬ 
mations  without  experiencing  issues  with  mesh  quality.  The  computations  presented 
so  far  had  no  body  motion.  We  consider  several  test  cases  to  observe  the  effect  of 
prescribed  body  motion  on  the  aerodynamic  forces  generated.  First  we  study  if  the 
small  amount  of  body  motion  seen  dnring  the  locnst  flight  has  an  effect  on  the  overall 
lift  and  thrust  production.  In  addition,  we  study  the  cases  of  larger  and  more  general 
prescribed  body  motion.  Specihcally,  we  consider  three  cases:  rolling,  pitching,  and 
combined  rolling  and  pitching. 
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8.1  Locust  body  motion 

We  prescribe  for  the  MAV  the  body  motion  that  the  locust  had  in  the  flight  from 
which  we  extracted  the  wing  motion  and  deformation  data. 

8.1.1  Setup  and  computations 

In  [3,  1,  82],  the  body  motion  of  the  locust  was  removed  from  the  motion  of  the  wing 
tracking  points.  Here  we  include  that  body  motion  separately,  which  consists  of  a 
translational-displacement  vector  and  three  rotations.  The  setup  for  the  representa¬ 
tion  of  the  wing  motion  is  identical  to  that  in  the  base  case.  In  the  temporal  repre¬ 
sentation  of  the  body  motion,  we  use  cubic  NURBS  basis  functions  for  the  position 
vector  and  the  three  angles.  With  this,  we  have  separate  temporal  representations  for 
the  body  and  wing  motions.  To  make  the  entire  motion  periodic,  the  body  position 
vector  and  three  angles  are  adjusted  just  like  how  the  wing  motion  was  made  periodic 
in  [1].  The  body  position  and  orientation  cycle  is  split  into  the  same  temporal  patches 
that  the  flapping  cycle  was  split. 

The  wing  flapping  is  handled,  as  before,  with  the  mesh  motion  inside  the  cylin¬ 
drical  rehnement  region.  The  body  motion  is  handled  by  moving  the  cylindrical 
rehnement  region,  which  in  turn  drives  the  motion  of  the  mesh  between  the  cylindri¬ 
cal  region  and  the  external  boundaries.  For  the  base  case,  there  is  no  mesh  motion 
outside  the  cylindrical  rehnement  region,  so  that  part  of  the  mesh  does  not  change. 
For  the  cases  with  body  motion,  because  the  locust  body  motion  is  fairly  small,  that 
mesh  changes  (i.e.  moves),  but  its  connectivity  does  not.  As  before,  we  subdivide  the 
prescribed  motion  into  smaller  steps  and  use  the  same  method  to  deform  the  outer  re¬ 
gion  as  we  did  for  the  mesh  motion  around  the  wings.  We  perform  the  huid  mechanics 
and  mesh  moving  computations  with  the  same  settings  as  the  base  case.  Figures  8.1 
and  8.2  show,  in  a  comparative  setting,  the  happing  and  body  motions  for  the  base 
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and  locust-body  cases.  Figures  8.3  and  8.4  show  the  vertical  and  longitudinal  wing 
tip  positions  and  velocities  for  the  base  and  locust-body  cases. 

8.1.2  Results 

Figure  8.5  shows  the  total  lift  and  thrust  for  the  base  and  locust-body  cases.  On  the 
downstroke,  the  lift  for  the  locust-body  case  is  somewhat  smaller  than  that  of  the 
base  case.  Also,  the  magnitude  of  the  thrust  is  smaller.  This  is  most  likely  due  to  the 
fact  that  on  the  downstroke,  the  body  motion  causes  the  body,  and  also  the  wings,  to 
move  up  and  back  thus  decreasing  the  actual  wind  velocity  around  the  wings.  This 
can  also  be  seen  during  the  upstroke,  with  the  magnitude  of  the  lift  and  thrust  being 
smaller  than  those  in  the  base  case.  Any  body  motion  opposite  the  motion  of  the 
wings  leads  to  a  decreased  relative  velocity  and  thus  decreased  force  magnitude.  The 
body  motion  for  this  case  is,  for  the  most  part,  opposite  to  that  of  the  wing  motion, 
and,  therefore,  leads  to  magnitude  loss  for  both  lift  and  thrust.  The  average  lift  and 
thrust  for  the  locust-body  case  is  6.04  mN  and  —1.27  mN,  respectfully,  compared 
to  6.50  mN  and  —1.31  mN  for  the  base  case. 


8.2  Prescribed  roll 

8.2.1  Setup  and  computations 

We  model  a  dynamic  roll  where  the  MAV  starts  at  18.5°,  rotates  to  45°,  then  down 
to  0°,  and  back  to  18.5°.  Such  a  large  range  of  motion  would  have  required  far 
more  frequent  remeshing  if  the  body  motion  drove  the  mesh  motion  directly.  Using 
the  external  mesh  to  take  up  the  motion  makes  the  mesh  moving  process  more  effi¬ 
cient.  The  large  range  of  rotation  requires  that  we  use  a  new  external  mesh  for  each 
patch.  However,  this  does  not  increase  the  remeshing  cost  that  much  (the  external 
mesh  has  roughly  30%  of  the  nodes).  The  HW  deformation  is  computed  with  the 
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Figure  8.1:  Flapping  and  body  motions  for  the  base  (grey)  and  locust-body  (blue) 
cases  at  the  first  four  of  eight  equally-spaced  instants  during  the  flapping  cycle.  Top 
(left)  and  side  (right)  views. 
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Figure  8.2:  Flapping  and  body  motions  for  the  base  (grey)  and  locust-body  (blue) 
cases  at  the  last  four  of  eight  equally-spaced  instants  during  the  flapping  cycle.  Top 
(left)  and  side  (right)  views. 
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Figure  8.3:  Vertical  tip  position  and  velocity  for  the  right  wings  for  the  base  and 
locust-body  cases. 
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Figure  8.4:  Longitudinal  tip  position  and  velocity  for  the  right  wings  for  the  base  and 
locust-body  cases. 
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Figure  8.5:  Total  lift  (top)  and  thrust  (bottom)  for  the  base  and  locust-body  cases. 
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Kirchhoff-Love  shell  model.  The  computational  parameters  for  the  fluid  and  struc¬ 
tural  mechanics  and  mesh  motion  are  the  same  as  in  Chapter  6.  Figure  8.6  shows  the 
body  and  wing  motions  used  in  the  fluid  mechanics  computation  of  the  hrst  SCFSl 
iteration.  We  perform  three  SCFSl  iterations  and  present  the  results  from  that. 

8.2.2  Results 

Figure  8.7  shows  the  total  lift  and  thrust  for  the  hxed-body  and  prescribed-roll  cases. 
The  rolling  motion  changes  the  aerodynamic  forces  generated  signihcantly.  The  pre¬ 
scribed  roll  results  in  a  noticeable  decrease  in  the  lift,  but  at  the  same  time  the  thrust 
production  is  improved  during  the  hrst  part  of  the  happing  cycle  (HW  downstroke). 
This  suggests  that  rolling  maneuvers  can  be  useful  in  some  cases  when  increased 
thrust  is  needed  and  lift  is  not  as  important. 

8.2.3  Prescribed  pitch 

Setup  and  computations 

The  setup  for  the  prescribed-pitch  case  is  very  similar  to  the  setup  for  the  prescribed- 
roll  case.  Instead  of  varying  the  roll  angle,  we  are  changing  the  body  angle,  which 
dynamically  alters  the  angle  of  attack  during  the  hight.  We  start  at  a  body  angle 
of  11.7°,  pitch  the  MAV  up  to  27.3°,  down  to  0°,  and  back  to  11.7°.  Again,  because 
the  rehnement  cylinder  has  a  large  motion,  a  new  external  mesh  is  used  for  each 
patch.  The  HW  deformation  is  compnted  with  the  Kirchhoff-Love  shell  model.  The 
compntational  parameters  for  the  flnid  and  strnctnral  mechanics  and  mesh  motion 
are  the  same  as  in  Chapter  6.  Fignre  8.8  shows  the  body  and  wing  motions  nsed 
in  the  flnid  mechanics  compntation  of  the  hrst  SCFSl  iteration.  We  perform  three 
SCFSl  iterations  and  present  the  results  from  that. 
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Figure  8.6:  Prescribed-roll  case  with  the  HW  deformation  extracted  from  the  struc¬ 
tural  mechanics  computation,  seen  at  eight  equally-spaced  instants  during  the  third 
flapping  cycle. 
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Figure  8.7:  Total  lift  (top)  and  thrust  (bottom)  for  the  hxed-body  and  prescribed-roll 


cases. 
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Figure  8.8:  Prescribed-pitch  case  with  the  HW  deformation  extracted  from  the  struc¬ 
tural  mechanics  computation,  seen  at  eight  equally-spaced  instants  during  the  third 
flapping  cycle.) 
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Results 


Figure  8.9  shows  the  total  lift  and  thrust  for  the  hxed-body  and  prescribed-pitch 
cases.  We  see  an  increase  in  both  the  lift  and  thrust  in  the  initial  part  of  the  HW 


Figure  8.9:  Total  lift  (top)  and  thrust  (bottom)  for  the  fixed-body  and  prescribed- 
pitch  cases. 


downstroke,  but  the  overall  lift  and  thrust  both  decrease.  During  that  hrst  part  of 
the  cycle  the  pitching  motion  is  upwards  (increases  the  angle  of  attack).  By  changing 
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the  direction  of  the  pitching  motion  with  respect  to  the  HW  motion,  we  can  influence 
the  total  lift  and  thrust  production  of  the  MAV. 

8.2.4  Prescribed  roll  and  pitch 

Setup  and  computations 

The  prescribed  body  motion  in  this  case  is  a  synchronized  combination  of  the  pre¬ 
scribed  roll  and  pitch  in  the  two  earlier  cases.  Again,  a  new  external  mesh  is  used  for 
each  patch.  The  HW  deformation  is  computed  with  the  Kirchhoff-Love  shell  model. 
The  computational  parameters  for  the  fluid  and  structural  mechanics  and  mesh  mo¬ 
tion  are  the  same  as  in  Chapter  6.  We  perform  three  SCFSI  iterations  and  present 
the  results  from  that. 

Results 

Figure  8.10  shows  the  total  lift  and  thrust  for  the  fixed-body  and  prescribed  roll  and 
pitch  cases.  At  the  beginning  of  the  flapping  cycle  there  is  an  increase  in  the  lift,  but 
the  overall  lift  production  is  decreased.  As  before,  the  average  lift  and  thrust  both 
decrease.  Figures  8.11  and  8.12  show  the  pressure  on  the  wing  surfaces  (relative  to 
the  free-stream  pressure)  at  eight  equally-spaced  instants  during  the  third  flapping 
cycle.  As  can  be  seen,  the  pressure  distribution  is  quite  asymmetrical  compared  to 


the  base  case. 


Force  (mN)  Force  (mN) 


102 


84 


t/T 


Figure  8.10:  Total  lift  (top)  and  thrust  (bottom)  for  the  hxed-body  and  prescribed 
roll  and  pitch  cases. 
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Figure  8.11:  Surface  pressure  in  Pa  (relative  to  the  free-stream  pressure)  for  the 
prescribed  roll  and  pitch  case,  seen  at  the  hrst  four  of  eight  equally-spaced  instants 
during  the  third  flapping  cycle  (top  view  on  left,  bottom  view  on  right). 


104 


86 


^ 

-14  0  14 

Figure  8.12:  Surface  pressure  in  Pa  (relative  to  the  free-stream  pressure)  for  the 
prescribed  roll  and  pitch  case,  seen  at  the  last  four  of  eight  equally-spaced  instants 
during  the  third  flapping  cycle  (top  view  on  left,  bottom  view  on  right). 
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Chapter  9 


Flapping-pattern  studies 


It  is  interesting  to  study  the  effect  of  changing  the  synchronization  between  the  FW 
and  HW  and  introduce  asymmetry  between  the  motions  of  the  left  and  right  wings. 
In  this  chapter  we  present  such  studies  from  [2],  In  the  first  study,  we  alter  the 
synchronization  between  the  FW  and  HW  by  advancing  or  delaying  the  FW  motion 
relative  to  what  it  is  in  the  original  locust  data.  This  allows  us  to  seek  an  optimum 
point  for  the  lift  and  thrust  production  and  evaluate  the  possibility  of  using  the 
variation  in  the  synchronization  to  alter  the  flight  performance.  In  the  second  study, 
we  introduce  asymmetry  (across  the  sagittal  plane)  to  the  flapping.  From  a  practical 
standpoint,  it  is  unlikely  to  have  asymmetric  flapping  in  the  MAV  design,  at  least  in 
the  current  ones.  Still,  it  would  be  valuable  to  have  the  ability  to  study  the  effect 
of  such  asymmetry  for  the  purpose  of  evaluating  more  complex  MAV  designs.  We 
created  a  test  case  where  the  left  and  right  wings  are  flapping  still  with  the  same 
frequency,  but  one  side  is  delayed  with  respect  to  the  other. 
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9.1  Setup  and  computations 

9.1.1  FW  advanced  or  delayed 

For  the  case  where  the  FW  is  advanced  (“FWA”)  or  delayed  (“FWD”),  the  setup 
is  very  similar  to  the  case  with  the  original  locust  data  (“FWO”),  other  than  just 
advancing  or  delaying  the  FW  by  9.1%  of  the  flapping  period.  The  FW  and  HW 
come  very  close  to  each  other  during  flapping,  so  it  is  not  possible  to  change  the  syn¬ 
chronization  that  much  because  the  wings  would  collide  unless  the  distance  between 
the  wings  along  the  body  is  increased.  As  this  would  cause  a  variety  of  other  changes 
and  make  the  comparisons  less  consistent,  we  use  a  value  close  to  the  maximum  pos¬ 
sible  advancing  and  delaying  that  does  not  require  increasing  the  distance  between 
the  wings.  Figure  9.1  shows  the  vertical  tip  position  of  the  FW  for  the  FWA  and 
FWD  cases  compared  to  the  FWO  case.  The  HW  deformation  is  computed  with  the 


Figure  9.1:  Vertical  tip  position  of  the  FW  for  the  FWA  and  FWD  cases  compared 
to  the  FWO  case. 


Kirchhoff-Love  shell  model.  The  computational  parameters  for  the  fluid  and  struc¬ 
tural  mechanics  and  mesh  motion  are  the  same  as  in  Chapter  6.  The  mesh  motion 
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strategy  in  the  refinement  cylinder  is  the  same  as  it  was  in  the  FWO  case.  We  perform 
three  SCFSI  iterations  and  present  the  results  from  that. 

9.1.2  Asymmetric  flapping 

In  this  case  we  are  delaying  the  left  wings  by  54%  of  the  flapping  period.  Figure  9.2 
shows  the  vertical  tip  positions  of  the  FW  and  HW  for  both  left  and  right  sides. 
The  HW  deformation  is  computed  with  the  Kirchhoff-Love  shell  model.  The  com¬ 
putational  parameters  for  the  fluid  and  structural  mechanics  and  mesh  motion  are 
the  same  as  in  Chapter  6.  The  mesh  motion  strategy  in  the  refinement  cylinder  is 
the  same  as  it  was  in  the  symmetric-flapping  case,  except  for  the  mesh  asymmetry 
between  the  two  sides  across  the  sagittal  plane.  Figure  9.3  shows  the  body  and  wing 
motions  used  in  the  fluid  mechanics  computation  of  the  first  SCFSI  iteration.  We 
perform  three  SCFSI  iterations  and  present  the  results  from  that. 

9.2  Results 

9.2.1  FW  advanced  or  delayed 

Figures  9. 4-9. 9  show  the  lift  and  thrust  for  the  FWA  and  FWD  cases  compared  to 
the  FWO  case.  The  average  lift  is  1.43  mN,  1.40  mN  and  1.39  mN  for  the  FWO, 
FWA  and  FWD  cases.  The  peak  lift  occurs  at  different  times.  The  average  thrust 
is  —0.56  mN,  —0.50  mN  and  —0.52  mN  for  the  FWO,  FWA  and  FWD  cases.  The 
original  synchronization  has  the  most  lift,  while  both  the  FWA  and  FWD  cases  show 


an  increase  in  thrust. 
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Figure  9.2:  Vertical  tip  positions  of  the  FW  and  HW  for  the  left  and  right  sides  in 
the  asymmetric-flapping  case. 
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Figure  9.3:  Asymmetric-flapping  case  with  the  HW  deformation  extracted  from  the 
structural  mechanics  computation,  seen  at  eight  equally-spaced  instants  during  the 
third  flapping  cycle. 
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Figure  9.4:  Total  lift  (top)  and  thrust  (bottom)  for  the  FWA  and  FWO  cases. 


Force  (mN)  Force  (mN) 


111 


93 


0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 

t/T 

Figure  9.5:  Total  lift  (top)  and  thrust  (bottom)  for  the  FWD  and  FWO  cases. 
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Figure  9.6:  Right  FW  lift  (top)  and  thrust  (bottom)  for  the  FWA  and  FWO  cases. 
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Figure  9.7:  Right  HW  lift  (top)  and  thrust  (bottom)  for  the  FWA  and  FWO  cases. 
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Figure  9.8:  Right  FW  lift  (top)  and  thrust  (bottom)  for  the  FWD  and  FWO  cases. 
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Figure  9.9:  Right  HW  lift  (top)  and  thrust  (bottom)  for  the  FWD  and  FWO  cases. 
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9.2.2  Asymmetric  flapping 

Figure  9.10  shows  the  total  lift  and  thrust  for  the  symmetric  and  asymmetric  flapping. 
Figures  9.11  and  9.12  show,  for  the  right  FW  and  HW,  the  lift  and  thrust  for  the 
symmetric  and  asymmetric  flapping.  Figures  9.13  and  9.14  do  that  for  the  left  FW 
and  HW. 
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Figure  9.10:  Total  lift  (top)  and  thrust  (bottom)  for  the  symmetric  and  asymmetric 
flapping. 
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Figure  9.11:  Right  FW  lift  (top)  and  thrust  (bottom)  for  the  symmetric  and  asym¬ 
metric  flapping. 
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Figure  9.12:  Right  HW  lift  (top)  and  thrust  (bottom)  for  the  symmetric  and  asym¬ 
metric  flapping. 
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Figure  9.13:  Left  FW  lift  (top)  and  thrust  (bottom)  for  the  symmetric  and  asymmetric 
flapping. 
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Figure  9.14:  Left  HW  lift  (top)  and  thrust  (bottom)  for  the  symmetric  and  asym¬ 
metric  flapping. 
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As  expected,  for  the  right  wings  the  lift  and  thrust  for  the  symmetric  and  asym¬ 
metric  flapping  are  very  similar.  For  the  left  wings  they  are  quite  different.  The  left 
HW  lift  and  thrust  are  basically  shifted  in  time  in  a  fashion  that  reflects  the  flapping 
shift  in  time.  The  left  FW,  on  the  other  hand,  has  very  different  lift  and  thrust 
prohles.  The  average  left  FW  lift  and  thrust  for  the  symmetric  flapping  are  0.63  mN 
and  0.01  mN.  For  the  asymmetric  flapping,  the  values  are  0.64  mN  and  —0.08  mN. 
The  total  lift  and  thrust  for  the  symmetric  flapping  are  1.43  mN  and  —0.56  mN.  For 
the  asymmetric  flapping,  the  values  are  1.40  mN  and  —0.68  mN. 
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Chapter  10 


Conclusions 


We  repeat  here  the  concluding  remarks  from  [2],  We  have  focused  on  a  SCFSI  anal¬ 
ysis  of  flapping-wing  aerodynamics  of  an  MAV.  The  wing  motion  and  deformation 
data,  whether  prescribed  fully  or  partially,  was  from  an  actual  locust,  extracted  from 
high-speed,  multi-camera  video  recordings  of  the  locust  in  a  wind  tunnel.  The  core 
computational  FSI  technology  is  based  on  the  DSD/SST  formulation,  specihcally, 
DSD/SST-VMST  (also  called  ST-VMS),  which  is  the  version  of  the  DSD/SST  for¬ 
mulation  derived  in  conjunction  with  the  VMS  method.  This  was  supplemented  with 
special  ST  techniques  targeting  flapping-wing  aerodynamics,  and  the  common  feature 
for  much  of  these  special  ST  techniques  is  using  NURBS  basis  functions  in  temporal 
representation.  That  includes  temporal  representation  of  the  wing  and  mesh  motion. 
Temporal  representation  with  NURBS  is  also  a  key  feature  of  the  remeshing  process. 
The  structural  mechanics  computations  were  mostly  based  on  the  Kirchhoff-Love 
shell  model.  We  also  carried  out,  for  comparison  purpose,  some  test  computations 
with  the  membrane  model.  The  sequential-coupling  technique,  which  is  applicable 
to  some  classes  of  FSI  problems,  especially  those  with  temporally-periodic  behavior, 
performed  well  in  FSI  computations  of  the  flapping-wing  aerodynamics  we  consid¬ 
ered  here.  In  addition  to  the  straight-flight  case,  we  analyzed  cases  where  the  MAV 
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body  had  rolling,  pitching,  or  rolling  and  pitching  motion,  where  we  assumed  that  the 
MAV  was  somehow  able  to  perform  such  maneuvers.  We  showed  how  these  maneuvers 
influence  the  lift  and  thrust. 
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